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The orbital degree of freedom plays a fundamental role in understanding the uncon¬ 
ventional properties in solid state materials. Experimental progress in quantum atomic 
gases has demonstrated that high orbitals in optical lattices can be used to construct 
quantum emulators of exotic models beyond natural crystals, where novel many-body 
states such as complex Bose-Einstein condensates and topological semimetals emerge. 
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mary of exotic orbital models and resulting many-body phases is provided. Experimental 
consequences of the novel phases are also discussed. 
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I. INTRODUCTION 


Optical lattices play a central role in studying strongly interacting many-body physics with ultracold atoms (Bloch 


et ai 2008 Dutta et ai 2015a Lewenstein et ai 2012). Because of their unprecedented controllability, atomic 
gases confined in optical lattices enable quantum simulation of various lattice Hamiltonians, e.g., Bose and Fermi 
Hubbard models, where different aspects have been intensively investigated. With single-species of bosons, e.g., ^^Rb, 
a quantum Mott-to-superfiuid transition has been observed. Multi-component lattice models have been reached with 
atomic internal degrees of freedom. The SU(2) spinfull Fermi-Hubbard simulator has been carried out by using 
hyperfine states of ^Li or One theme along this direction is to emulate complex correlated phenomena of strongly 
interacting electrons. Such multi-component quantum simulators with atomic internal degrees of freedom have been 
very successful in simulating Hamiltonians with high symmetries. 

For electrons, one important ingredient besides spin is the orbital degrees of freedom, which arises in a variety of 
condensed matter systems (Tokura and Nagaosa 2QQQ| ). In solid state materials, orbitals originate from electron clouds 
surrounding the ions in the crystal. With tunnelings, these orbitals form Bloch bands. Orbitals are Wannier states 
corresponding to different bands. Degenerate orbitals (or bands) could emerge in presence of point group symmetries, 
but the symmetry for orbitals is much lower than for spins. Understanding such orbitals degree of freedom is crucial 
to obtain a simple and yet powerful model that captures the essence of many complicated materials, such as transition 
metal oxides, pnictides, etc. A task of this kind however remains outstanding, much due to the complexity of multiple 
types of degrees of freedoms coupled together, including orbital, charge, spin, and crystal field. The intricate coupling 
makes it an expensive challenge, both analytically and numerically, to understand orbital physics alone first and to 
attempt to compare with any electronic solid state materials in experiments. 

Given one important application of optical lattices is to simulate complex phenomena of electrons, it is rather 
essential to find ways to emulate electron orbitals with atoms. Actually with optical lattices, the ionic crystal trapping 
electrons is replaced by an artificial crystal of light, created by standing waves of laser beams. The Wannier orbitals 
in the lattice naturally mimic properties of that in ionic crystals. Due to the intrinsic spatial nature, orbital degree of 
freedom in both of these ionic and light crystals respect space point group symmetries rather than internal continuous 
group symmetries, which defines its uniqueness. Such symmetry properties of orbitals make them fundamentally 
difficult to be simulated with internal atomic degrees of freedom such as hyperfine spins. On this regard, the orbital 
states of an atom in an optical lattice provide a natural avenue to emulating the electronic orbital related physics. 

Exploration of orbital physics in optical lattices is certainly not restricted to quantum simulations of electrons in 
solids. For example, orbital bosons are able to bring to the study of quantum matter some really novel concepts that 
have no prior analogue in systems of (fermionic) electrons. Moreover, bosons (e.g., ^^Rb atoms) are more widely used 
in optical lattice experiments. In the first experimental demonstration of many-body orbital physics, bosons were 
loaded into the p-bands of an optical lattice, for which earlier theoretical studies had predicted interesting phenomena 


such as time-reversal symmetry breaking and spontaneous angular momentum order (Isacsson and Girvin 2005 


Kuklov 2006 Liu and Wu 2006). 


Strong interactions which are achievable in optical lattice experiments also lead to interesting orbital physics. 
Firstly, with strongly repulsive bosons loaded into higher orbital bands, they would form a Mott state with orbital 
degree of freedom. Orbital ordering in such a Mott state is drastically different from spin ordering in Mott states. For 
Mott states formed by spinor bosons (assuming no spin-orbital coupling), the super-exchange Hamiltonian typically 
has high symmetries. The orbital super-exchange Hamiltonian is generally more complicated and at the same time 
promises richer physics. Secondly, for strongly interacting atoms in a lattice (e.g., lattice bosons in the Mott regime, 
or a Feshbach resonant Fermi gas (Ghin et al 2010) in a lattice), even without deliberately loading atoms into higher 
orbital bands, population of those bands is unavoidable due to interaction effects. This is because local interactions 


would mix all different orbitals. Recent works (Soltan-Panahi et al. 2012 Zhou et al. 2011b) have shown that the 


interaction-induced high-band population could give rise to significant physical effects, such as condensation of boson 
pairs, and exotic symmetry breaking orders. It is therefore essential rather than an option to account for orbital 
physics in modeling strong interaction effects in optical lattices. 

Research of fermions in higher orbitals adds a remarkably distinct venue. Theoretical studies have also found 
quantum phases with angular momentum ordering that spontaneously break time-reversal symmetry. For fermions, 
this symmetry breaking leads to even more dramatic effects than the bosonic counterpart. Gonsidering the angular 
momentum order and mixing of orbitals with opposite parities (like s and p, or p and d orbitals), the fermionic 
atoms experience effective gauge fields, which then gives rise to topological phenomena, like quantum Hall, topolog¬ 
ical insulator, or certain topologically protected gapless phases. This route of engineering topological matter offers 


one way different from the Raman-induced synthetic gauge fields (Dalibard et al. 2011) or the artificial spin-orbit 


couplings (Galitski and Spielman 2013 Zhai 2015). It has fundamentally distinct properties and is advantageous 
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in certain aspects. For example, it does not involve complications of Raman couplings, and the resultant topological 
phases would have longer lifetime due to less heating effects. The finite temperature behaviors of the spontaneously 
generated gauge fields are also different from the the Raman-induced case. 

In this review, we start by describing basics of modeling orbitals in optical lattices. Then by using particular 
examples, we present a selection of many-body aspects of orbital physics that we find most interesting and novel, 
as sketched above. Along with developing theoretical concepts and models pedagogically, we review the recent 
experimental developments and the current status in this field, and outline several future directions. 


II. HIGH ORBITALS AND BAND STRUCTURES IN OPTICAL LATTICES 


Previous studies in optical lattices largely focused on atoms trapped in the lowest band and the resultant single¬ 
band Hubbard model, where correlated effects of bosons, e.g., the Mott-superfluid transition, have been intensively 


investigated (Bloch et ai 2008 Dutta et aL, 2015a Lewenstein et al. 2012). In this section we present the procedure 


to construct tight binding models involving high orbital degrees of freedom, which is one essential step to study 
correlation effects in interacting atoms in lattices. To demonstrate the validity condition of tight binding models, we 
also show the exact results from plane-wave expansion for the tunneling amplitudes, band structures and Wannier 
functions of higher bands. A two dimensional square lattice is assumed in this section. 


A. Harmonic approximation and tight binding models 


In the tight binding regime, an optical lattice can be treated as individual harmonic oscillators, which are coupled by 
quantum tunnelings. On each harmonic oscillator centered at a lattice site labeled by its position R, we have discrete 
energy levels with orbital wavefunctions 0 q,(x —R). Associated with the localized orbital wavefunctions, we can define 
the lattice operators 6 q,(R). To do this, it has to be enforced that the orbital wavefunctions are orthonormal. The 
simple eigen wavefunctions of harmonic oscillators do not satisfy orthonormal condition, for the reason that there are 
overlaps between orbital wavefunctions on neighboring sites. 

The procedure to construct the orthogonal basis from the localized harmonic oscillator wavefunctions is the following. 
We start with the harmonic oscillator wavefunctions (pa (x — R) localized on site R. These wavefunctions are already 
approximately orthogonal, i.e.. 


/ 


dx(p*^{x - R)(/)a/(x - R') = + e„R,a'R', 


where eaii,a'ii' are small numbers. By definition we know that [e] is a traceless Hermitian matrix. Then we improve 
this basis by introducing 


^(x) = (f>a{x) - ^ y] ea'R',aR</>a'(x - R') 


( 1 ) 


a'R' 


After that <pa (x — R) is renormalized as 


0a(x) ^ (^c«(x)/-y/ J (ix'|0„(x')p. 

The improved wavefunctions satisfy a better approximate orthogonal condition 

J c^x^«(x - R)^a'(x - R') = <5aa'<5RR/ + O(e^). 

The above procedure can be iterated N times to get the orthonormal basis to the precision of 
Once we have the orthonormal basis, the tunnelings between R and R' are calculated as 

t„c«'(R-R') = y'rfxC(x-R)R(x)0„'(x-R'), (2) 

where i^(x) is the Hamiltonian in the first quantization form H = + ^(^)- The lattice model Hamiltonian 

including tunnelings is given by 

H= W'(R-R')&l(R)^c«'(R')- 

aa', RR' 


(3) 
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Without truncating the basis, the Hamiltonian is exact, from which the band structure can be calculated. If we only 
keep the lowest harmonic wave functions, this lattice Hamiltonian gives qualitatively correct band structures for deep 
lattices. 

The procedure described above to construct orthogonal orbital wave functions is one essential step if one uses 
harmonic approximation. In principle, the constructed wave functions are not the same as the maximally localized 


Wannier functions (Kivelson 1982 Marzari a/. 2012 Uehlinger a/. 2013 Ganczarek et a/. 2014). The procedure 


to calculate such maximally localized Wannier functions is not as straightforward and is beyond the scope of this 


review. 


Multi-band Hubbard model — Considering interacting bosonic atoms loaded on excited bands, the physics will be 
described by a multi-band Hubbard model 


H=J2 - a')bi{a)ba{a') 

RR' 

T ^ ^ b^QlQ2Q3Q4^Ql (^)^Q2 (^) * 

R 


With weak interaction, the coupling constants U^ 


et al 1998 Li et al. 2014a Liu and Wu 2006 Zhou et al. 2011b) 


QiQ 2 Q 3 Q :4 


can be estimated at tree-level as (Dutta ei aL 2015a Jaksch 


U, 


aia2Ct3Ct4 


2m 


/ 


(x)(/>* 2 (x)(/>a3 (x)(/>a4 (x), 


( 5 ) 


where is the s-wave scattering length, tunable with Feshbach Resonance techniques. With fermionic atoms, we 
have a similar Hubbard model with interactions between hyperfine states. 


B. Band structures 


In terms of field operators, the Hamiltonian of particles moving in optical lattices is 

H = I + l/(x)^ t/>(x), (6) 

where '0(x) is a field operator. It can be either bosonic or fermionic. Statistics is irrelevant here to determine 
single-particle band structures. We expand the operator '^(x) in the momentum basis 

^(x) = E ^ E «K(k)e*(K+k)-x^ (7) 

where K labels the reciprocal lattice vectors, k the lattice momentum, and Ns the number of lattice sites. Here and 
henceforth, the lattice constant is set as the length unit. Optical lattice potentials H(x), unlike the potentials in 
electronic materials, can typically be written as superpositions of just a few plane waves, i.e.. 


V{x) = 

K 

For example, the potential of a square lattice created by laser is 

'F(x) = —Vq [sm‘^{kx) + sin^(%)] 

= -^ -k C.C.] -k const, 

where k is the wavevector of the laser beams. The Hamiltonian in momentum space reads as 

^ = E E ^k(Ki,K2)4^(k)dK,(k), (8) 

k Ki,K2 

with the matrix given by 

nk(Ki,K2)=^^l^^^NlsK^K,+v{Ki-K2). (9) 
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Diagonalizing this matrix, we get the band structure En(k) and the eigenvectors A^^k), with n the band index. The 
Hamiltonian in the eigen-basis reads 


ff = EE^"(k)5t(k)5„(k), 

n k 


( 10 ) 


with 6„(k) = A^^*(k)aK(k). 


VqIEr AtnniER AtnnniER W^nlER At^rinlER 


3 

-0.4441 

0.0449 

2.0074 

0.3308 

5 

-0.2631 

0.0136 

1.6912 

0.2914 

10 

-0.07673 

9.1E-4 

0.9741 

0.1051 

20 

-9.965E-3 

1.2E-5 

0.2411 

5.5E-3 


TABLE I Tunneling amplitudes in a two dimensional square lattice with potential E(x) = — Vb [sin^(/ca:) +sin^(/cy)]. Er is 

the one photon recoil energy tnn cind t^nn sire nearest neighbor and next nearest neighbor tunnelings for the lowest s 

band, and Einn nearest neighbor and next nearest neighbor tunnelings in the x direction for the px (first excited) band. 


The Wannier basis is given by 

= ( 11 ) 

Inversely we have 6n(k) = 

The Wannier wavefunctions of the Bloch bands are given by 


w. 


.(x-R) = E 


K 




The Hamiltonian can be rewritten in the Wannier basis as 

R = E - n')bl(R)bn(R'), 


RR' 


with 




( 12 ) 


(13) 


(14) 


Typical values of tunnelings (tunnelings refer to tunneling matrix elements here) for s and p bands are listed in Table 
In the definition of Wannier functions (Eq. (p!^), there are gauge degrees of freedom A^^(k) ^ e*^^^^^A^^(k). One 
has to make a smooth gauge choice to get localized Wannier functions (Marzari et al. |2Q12 ). Wannier functions of 
p-bands of a square lattice are shown in Fig. 


C. Heuristics to lifetime of high orbital atoms 


Here the lifetime of p-orbital condensate in a one dimensional (ID) lattice is discussed based upon Fermi’s Golden 
rule calculation. The resulting time scale is expected to apply to two dimensional (2D) square and three dimensional 
(3D) cubic lattices as well (Isacsson and Girvin 2QQ5[ ). The p-orbital condensate wavefunction is given as 


1 AT 


Vn\ 


(15) 


With interactions, two particles in the p-band can collide and one particle would decay to the lowest s-band and the 
other goes to the second excited d-band. This process is described by the following interaction term 


HZ^ = U/N, E {blih)bl,{k2)bpih)bp{h) + h.c.}. 

ki-\-k2+ks-\-k4=0 
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FIG. 1 Lowest two bands of a square lattice. The potential we choose here is V (x) = —Vo [sin^(/cx) + sin^ (ky )^, with VqIEr = 4 
(Er is the one photon recoil energy), (a) shows band structures of px and s bands, whose Wannier functions are respectively 
shown in (b) and (c). 


The final state after the collision is 

= &J(fci)6]((fc2)J-^===|vac). 

The transition probability from second order perturbation theory is 

p{k„k2-,t) ^ 

where Ae/c^/c 2 difference of kinetic energy between |^) and /ci, A: 2 ). The loss rate from the p-band is obtained 


wt = E 


t 

kik2 

k\k2 


hN, 


-N{N-1) 


1 


1 


n -1 


p{e,iK)) pie^i-K)) 


where p(e) is the density of states and K is determined by 

es{K)^ed{-K) = 2ep{Q), 

which in general has two solutions when the band gap between s and p matches that between p and d. 
The loss rate per site is 


A7t{uUY 


1 


1 


p{es{K)) p{e,{-K))\ ’ 


-1 


(16) 


(17) 
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with V the filling factor. The lifetime Xjw is typically short for cubic or square lattices, where the condition of 


Eq. (16) may be satisfied. It was suggested that anharmonicity (Miiller et al 2007) present in the actual optical 
lattice potential should help suppress the decay. Nonetheless, the lifetime can be significantly improved by using 


double-well lattice potentials to mismatch the band gaps as first discussed in Ref. (Stojanovic et al. , 2008) and further 


confirmed in the experiments (Wirth et al 2011). 


III. MANY-BODY PHASES AND TRANSITIONS 


Orbital degrees of freedom play an important role in understanding many complex phases in solid state materials. 
For example, high temperature superconductivity in the cuprates ( [Bednorz and Muller 1986) and pnictides (Kami- 


hara et al. 2006), chiral p-wave superconductivity proposed in Sr 2 Ru 04 (Luke et al. 


ated by strong correlation effects in a multi-orbital setting (Tokura and Nagaosa 2000). In optical lattices, recent 


1998), and Ferromagnetic 


superconductivities in oxide heterostructures such as LaAl Os/SrTiOa (|Ohtomo and Hwang 2004), are all nude- 


studies have shown that the interplay of high orbitals and interaction effects give rise to unconventional many-body 


phenomena (Lewenstein and Liu 2011). 


For bosons loaded into high-orbital bands of an optical lattice, an analogue of Hand’s rule coupling leads to a complex 


Bose-Einstein condensate with spontaneous angular momentum order (Isacsson and Girvin 2005 Kuklov 2006 Liu 


and Wu 2006 Wu 2009 Wu et al. 2006). The bosonic analogue of Hand’s rule basically states that repulsive 
contact interactions favor maximization of the local angular momentum. Different aspects of the unconventional 
condensate have been theoretically investigated, e.g., rotation effects (Umucahlar and Oktell |2008|), manifestations of 


lattice geometry and trapping potential (Cai and Wu 


transitions (Pietraszewicz et al. 2012 Pinheiro et al. 


2011 


2013 


Lim et al. , 2008 Pinheiro et al. 


2015 Stasyuk and Velychko 2011 2012). Experimentally, 


20121), and orbital phase 


this complex Bose-Einstein condensate has recen tly been demonstrated in a checkerboard optical lattice (Kock et al. 


2015 

Olschlager et al. 

2013 

Wirth et al. 

2011) 


2011). By considering strong interactions, this condensate state develops 


a quantum phase transition to a Mott state with very rich orbital ordering, which has been studied by mean field 


theories (Collin et al.\ 2010 Larson et al. 


et al. , 2013 Li et al. , 2012 Sowinski et al. 


2009 


Li et al. 2011b) and also by unbiased numerical methods (Hebert 


2013). Even without deliberately loading atoms into the higher bands, it 


has been shown high-band population can be stabilized by interaction effects (Soltan-Panahi et al. 2012 Zhou et al. 


2011b). 


For fermions, it has been shown that interaction effects combined with the band topology of p-orbitals lead to 
various exotic quantum phases. With p-orbital fermions in two dimensions, interactions cause generic instabilities 


towards quantum density wave orders (modulations in spin, charge or orbital 

density) (|Lu and Arrigoni 2009 Wu 

2008b Wu and Das Sarma 

20081 Wu et al. 

2006 

12012 Zhang et al. 2012| 

Zhao and Liu 2008), unconventional 

Gooper pairings (Gai et al. 

2011 [Hung et al. 

2011 

|Lee et al. 2010 Liu et al.\ 

|2015| 

Zhang et al. 

2011 

2010b), and 

novel quantum magnetism (Hauke et al, 2011 

Wang et al. 2008 Wu and Zhai| 

|2008| 

Zhang et al. 

2010a 

Zhou et al. 


2015 ) at low temperature. From quantum engineering perspectives, the elongated spatial nature of p-orbitals makes 
them ideal building blocks for fascinating topological states, e.g., topological semi-metal ( |Sun et (Z',|2012b ), quantum 


2016 2014 2010), and even fractional states (Sun et a/. 2010[ ). 


Hall phases (Wang and Gong 2010 Wu 2008a), topological insnlators/snperconductors ( |Li et d. 2013 Liu et al. 


In this section, we will review a selection of quantum many-body phases of p-orbital bosons and fermions. 


A. Orbital pip Bose-Einstein condensation 


1. Complex Pa; ipy Bose-Einstein condensation at finite momentum 


For bosons loaded on p-orbitals of a 2D square lattice (Wirth et al. 2011) in the tight-binding regime, the tunneling 
Hamiltonian is 


Htun — E {^11 [bU^)bx{r + ax)+x ^y] 

r 

- t± [bl{r)b^{r + ay) + x y] + h.c.} , (18) 

where and by are bosonic annihilation operators for and Py orbitals, respectively (Fig. |^. After a Fourier 
transformation, we get the energy spectra for the Px and Py bands. The dispersion for the Px band is 


ea,(k) = 2t|| cos{kx) — 2t±_ cos{ky). 


















































































































































































































































8 

The dispersion for the Py band is readily obtained with a lattice rotation (C 4 ). There are two degenerate minima— 
= (tt, 0) and = (0, tt) in the p-bands with the degeneracy protected by the C 4 symmetry. The ground state 
manifold of non-interacting p-orbital bosons is spanned by 


\N,.Ny) 


[4(QJ] ^ [^^(Q,)] 


JV„ 


I vac), 


(19) 


which has a large degeneracy that shall be lifted by interactions. 



FIG. 2 Illustration of the tight binding model of p-orbital bosons on a square lattice (Liu and Wu 2006). The longitudinal 
tunneling amplitude ty is in general far greater than the transverse tunneling t±. The “±” symbols indicate the sign of two 
lobes of p-orbital wave functions. 


The interaction terms of repulsive p-orbital bosons read (Isacsson and Girvin 2005 Li et al. 2011b Liu and Wu 


2006) 


^int = 


E 


-Ui [n,r(r)n,^(r) +nj^(r)nj^(r)] 


+2U-2,nx{v)ny{Y) 

+ \uz [bUr)bl{r)by{r)by{r) + h.c] |, 


( 20 ) 


where the density operators Approximating Wannier functions by localized harmonic wavefunctions, we 

have 


Ui = 3 U 2 =3U3 = U>0, 


( 21 ) 


from which the interaction can be rewritten as 




R ■- 


^(R) - -Ll{R) 


( 22 ) 


with n = blbjy and = ib^by + h.c. We thus expect that the angular momentum order is “universally” favorable 
in p-orbital Bose gases. 

It is however worth emphasizing here that the angular momentum ordering does not rely on the strict equality in 
Eq. ( 21 ) or the interaction form in Eq. ( 22 ). This becomes more clear with Ginzburg-Landau or effective field theories 

2013). Detailed studies taking into account unharmonic corrections and 


analysis (Li et al. 2014a 2012 


Liu et al. 


trapping potentials also confirm that the angular momentum order indeed exists in the regimes accessible to optical 


lattice experiments (Collin et al.\ 2010 Pietraszewicz et al. 2013 Pinheiro et al. 2012 Sowinski et al. 2013). 


To capture quantum/thermal fluctuations, two slowly varying bosonic fields are introduced as 


<A^(x) = E^a 


(Q^ + k)e^ 


ikx 
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where A is a momentum cut off. The effective Hamiltonian of the field theory of 0^(x) is 

H = J(fr [Ki {da;<pl{r)d:,(j):,{r) + x ^ y) 

+K 2 {dy4)l{r)dy<p^{r) + x y) 

-M +X^ y) 

+ 2^1 {4’l:4>x4’l:<Px + a; y) + 2g2<l)l.(l)x4'l<l>v 

+ \9i{4'l4>l<t>y<Py + h.c.) . 

In a superfluid state, we have {(j)y) = (fy. At mean field level, the energy of this state is 


(23) 


^ — 2^1 {\^x\^ + + ‘^g2\^xf\^yf + 2^3 + C.C.) . 

From Eq. ( [2T] ), we have gi = 3^2 = 3^3 > 0, and the relative phase between and Py is locked at ±7r/2, i.e. px = 
, where the “±” sign is spontaneously chosen. The superfluid state has a staggered angular momentum order 
{ — {LzCR)), which breaks time-reversal symmetry. Such a superfluid state is named transversely staggered 

orbital current (TSOC) superfluid. The phase configuration of this superfluid state and its momentum distribution 
are shown in Fig. 




2 4 


kx in units of n/a 


■High 

I 

■Low 


FIG. 3 Transverse staggered orbital current superfluid state (Li et al 2011b Liu and Wu 2006). (a) shows the phase 

configuration, from which one can infer the orbital current alternates from site to site, (b) shows the momentum distribution, 
which is confirmed in the experiments (Wirth et al 2011). 


An alternative way of looking at time-reversal symmetry breaking is to project interactions into the subspace 
spanned by \NxNy) in Eq. (19). In this subspace, the interaction reads 

(24) 


= + ^y)^N^N'jNyN^ + ^^NxNySjy^N'jNy 


+ 


Us 


2K 

-\-x y 

For the two orbital components to be miscible, we need 


\l^'xi^'x - U)^v{Ny - 1 ) N^+2^N'yNy-2 


2 U 2 - \Us\ < Ui 


( 25 ) 
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as analogous to spin miscible condition in spinor condensates (Pethick and Smith 2008). 
The angular momentum correlation is given by 

~{'^x){{'^y) T 1) ~ {'^y){{'^x) + 1) 


(26) 


with (...) the ground state expectation value. With Us > 0, to minimize the energy in Eq. (24), / {bl{Qx)^y{Qy)) + 


gets a negative value in the ground state and in the thermodynamical limit {{ux) ^ 1 , {riy) ^ 1 ), it approaches 
{—)2{nx){ny). The system thus has a long range correlation in angular momentum, i.e., 

const 7 ^ 0. The corresponding Ising order parameter is a staggered angular momentum Lz{H) = {—l)^^~^^yLz(R). 
When Us is negative, {{bl{Qx)by{Qy)) + h.c.) becomes positive, and the angular momentum order (Lz) vanishes 
and the system develops the other Ising orbital order Px -^Py with an order parameter (—(6|,(R)6^(R) + h.c.). 
From the above analysis, the transition at [/s = 0 is predicted to be first order (Fig. [^, although fluctuations may 
stabilize some intermediate state and the first order transition could be replaced by a sequence of double second order 
transitions. 


1 

'SO.5 


0 


-1 -0.5 0 0.5 1 

UUU2 

FIG. 4 Staggered angular momentum order. For positive t/ 3 , the staggered angular momentum order is finite and the condensate 
has a staggered px ± ipy (TSOC) order; while for the negative case, the staggered angular momentum order vanishes and the 
condensate has a px i Py order. In this plot, we assumed 2 U 2 — \U 3 \ < Ui such that the orbital mixed state has lower energy 
than Px or py state. 



2. Symmetry based effective field theory description 


The predicted TSOC superfluid state in the p-band tight binding model is also confirmed with effective field theory 
(EFT) treatment ( |Li et al. 2014a), which infers that the TSOC superfluid does not necessarily require a deep lattice. 
In the band structure calculation for the case of lattice rotation symmetry ( jWirth et al. 2011), dispersion of the 
relevant p-band Ep{k.) has two degenerate minima at = (tt, 0) and = (O^tt), around which low energy modes 
can be excited due to quantum or thermal fluctuations. This leads to a two-component EFT, where the fields are 
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introduced as 


.(X) = / 


d?(\ 

(27r)' 


rKQa+9)e" 


(27) 


with A a momentum cutoff and h{Qoc-\-q) annihilation operators for the Bloch modes near the band minima. The form 
of EFT is determined by considering lattice rotation and reflection symmetries, under which the fields transform 
as 


and 


4>x{x,y) 

_V 

-4>y{y,-x) 

. (t>v{x,y) _ 

7 

<Px{y,-x) 

4>x{x,y) 

_V 

1 

T 

_ 4>v{x,y) _ 

7 

(py{-x,y) 


(28) 

(29) 


respectively. The Hamiltonian density of the EFT consistent with these symmetries is 

/ ^2 g2 \ 

Heff = 4(x) ^ ^ M^) +x^y 

+ E 9a^a2\(t)a^?\4>a2? + 9s + h.c) , 


(30) 


with effective couplings i^n, and ^’s. This form of EFT is solely symmetry based, i.e., independent of microscopic 


details. For weakly interacting bosons, the coupling constants in Eq. (|30|) can be calculated from microscopic models 
(see Appendix 0- 


In the vicinity of thermal phase transitions of the superfluid phases, classical phase fluctuations are expected to 
dominate the universal physics, which allows us to ignore the sub dominant density fluctuations and to replace by 
with p the total density. In terms of phases the Hamiltonian density is rewritten as 

Heff = - K^{dye,f+x^ y) 




+ ^93P^ cos ( 2 ( 6 »a; - By)) . 


(31) 


Bearing in mind the periodic nature of the phases 0^^ a proper lattice regularization of this EFT leads to a coupled 
XY model. 


-f^phasiX [{^'^11 cos(A^ 6 >^(r)) - 2Jj_ cos(Ay 6 >^(r))} 

r 

+ {x^y}]-U y]sin2(6»a;(r) - 6 »y(r)), 


(32) 


where Aj^Q,(r) = ^Q,(r + aj) — ^Q,(r) with j = x, y. At zero temperature we have the TSOC superfluid where the phases 
are locked at 0x{^) = + 6 >o, 0y{Y) = r^ 7 r + 6 >o + 57 r/ 2 , with s = ± and Oq G [ 0 , 27r) spontaneously chosen. At finite 

temperature, the coupled XY model supports two types of topological defects. The first is a vortex in the phase 6 >o, 
which is a point defect with logarithmic energy cost. The second is a domain wall connecting two Ising domains with 
different s. Upon heating the TSOC superfluid, vortex proliferation should drive a Kosterlitz-Thouless transition and 
the domain wall fluctuations should drive an Ising transition. Monte Carlo study finds that the Kosterlitz-Thouless 
transition temperature is lower than the Ising transition ( Li et al\ |2QI4a ). 

From the effective field theory analysis, the p-orbital angular momentum order, or equivalently the phase 


locking, does not rely on the precise form of the interaction (Eq. (22)). The requirements are ^3 > 0 and two p 
orbitals being miscible. 


3. Population of higher bands by interaction 

Here, we will focus on condensation of weakly interacting bosons in a lattice potential. With weak interaction, the 
condensate is well described by Gross-Pitaevskii approach where the condensate wavefunction 0(x) is obtained by 
minimizing an energy functional 

£^GP = j {x.) f + ^(x) - M I </>(x) + 5 r|</>(x)|^, 


( 33 ) 
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With infinitesimal interaction, the condensate wavefunction resembles the lowest band Bloch wavefunction with lattice 
momentum k = 0, i.e., (/)(x) oc — R). This wavefunction preserves lattice translation and time-reversal 

symmetries meaning 0(x) = (/)(x + a) and 0 = 0*. From these preserved symmetries the generic form of the 
condensate wavefunction with weak interaction is 




(34) 


R 


provided that there are no first order transitions. The coefficients are real and the interaction induced high band 
condensate is at zero lattice momentum. In terms of A^, the energy Eqp reads as 


^GP = ^ (^n(k = 0) - /i) A^ + UooOo\ 


+4 ^ ^ b^OOOn^O^n + 10 E m^n 

+O{xl^o), 

with the interactions Uaia 2 a 3 a 4 introduced in Eq. (§. Minimizing this energy functional leads to 

2 _ IJ- - Eo{0) 

2f/oooo ’ 

2 ^^ 00071-^0 


^n >0 


En{0) - At' 


The ratio is readily given as 


Ao 


Uooon E'o(O) 

1^0000 .^n(0)—/i_ 


(35) 

(36) 

(37) 

(38) 


Physically, the high band condensate is due to competition of interaction energy and lattice potential energy—the 
interaction favors an extended condensate; while the potential energy favors a condensate with high density at the 
lattice minima. The high band population due to interaction effects is found in various settings (Alon et aL\ [2005 


Dutta et aL\ 20II||Guo et al., 2007||Hofer et al. 

poT^ 

1 Kantian et al.\ 2007| [Lacki et al.\ 20I3| 

Larson et al. 

2009 

Lhmann et al.\ [20121 Mering and Eleischhauer 

20111 

van Oosten et al.\ [2003| [Sakmann et al. 

[20II| [Zhou et al. 


In experiments, one can make a large fraction of high band condensate with double-well lattices, where the gap 


between lowest two bands is typically small and the fraction can be measured with band mapping techniques (Greiner 


et al. 2001). Experimental evidence of this phenomenon has recently been achieved (Soltan-Panahi et al 


2012 ). 


Eurthermore, on general argument, the double-well lattices can have the energy gap between the ground 5-band and 
the first excited p-bands significantly smaller than that between the first excited bands and higher bands (e.g., between 
p and d). This mechanism suppresses the decay by energy conservation law, making the first excited bands effectively 
metastable (Stojanovic et al. 2008). This will be further discussed in Section [Tv} 


Eor parity-symmetric lattices, the condensate wavefunction is parity even—0(x) = 0(—x), which implies Xuodd — 0, 
riodd referring to the parity odd bands with renodd(^) = ~'^nodd(~^)* mean field level, parity odd bands do 


not contribute and 

= 0. 

However they can form pair condensate orderings— (6^(r)6^/(r)) due to Gaussian 

fluctuations ( 

o 

ISl 

2011a|b 

). The mean field state is given by M) = exp (J(i^x0(x)'0'*'(x)) vac). The effective 

Hamiltonian of high banc 

modes at Gaussian level reads 


+ E [U00nmXlbl(k)bli-k) + h.C.]. 

nm 7 ^ 0 ,k 7^0 

Erom standard perturbation theory, the correction on the mean field state from high band fluctuations is 

UoOnm XlbUk)bl{-k)\M) 


E 

nm^^Ojk 


En{k) + -E'm(-k) - 2/2 


(39) 


( 40 ) 
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which mediates pairings 


(^n,r^m,r) — ^ ^ k)) 


/ 




U^OOnm^O 


{2 'itY £J„(k) + -E'TO(-k) - 2 //' 


In parity odd bands, bosons form pair condensate with {bnodd,r) — ^ 


^odd ^odd 


,r) Y 0- 


(41) 


4. Three dimensional p-orbital BEC and frustrated orbital ordering 


For bosons loaded on p-bands of a three dimensional cubic lattice, the tight binding Hamiltonian is (Liu and Wu 


2006) 


a — ^^[t\\^a/3 ^±(1 ^a/ 3 )] 


ra^ 

-Is 


=-) 


2 ^ t2 

'^r - o-^r 
o 


(42) 


where n and L are boson density and angular momentum operators and L^r = 

Without interaction, there are three degenerate p-bands and the energy minima are at = ( 7 r, 0 , 0 ), = ( 0 , 7 r, 0 ) 
and = (0,0, tt). The degenerate single-particle states are jQ^,) = 6 j^(QQ,)|vac). Thus any condensate wavefunction 
of a linear superposition of jQ^,) (Cai et ai | 2012 b| |, 




has the same single-particle energy. Here c = {cx^Cy^Cz) is a complex vector normalized to 1, i.e., |^ = 1. This 
complex vector could be parametrized as ( Liu and Wn[ [2006 ) 


Cy 


8 

_1 

1 

1 _ 


1 

0 


(43) 


^a=x,y,z the generators of SO(3) orbital rotation in the following matrix representation: 


with Ta 

Although the SO(3) orbital rotation is not a symmetry of the total Hamiltonian, it keeps the interaction term 
invariant because Ur and are both SO(3) scalars. With a condensate at the single-particle state |^, the mean field 
interaction energy is readily given as (Liu and Wu 2006) 


Eint = ^UNsnl 


1 - 3 sin2(2x) 


(44) 


with no the boson occupation number per site. The interaction energy is minimized at y = Similar to the two 
dimensional case, the time-reversal symmetry is spontaneously broken in the p-band condensate. The ground state 
manifold is U{1) x Z 2 x SO{3) at mean field level. The Z 2 x U{1) degeneracy remains due to the symmetries of the 
Hamiltonian, whereas the SO(3) degeneracy is an artifact of the mean field theory and such a degeneracy is lifted by 
fluctuations through an “order by disorder” mechanism. Ref. ( |Cai et oL 2012 b) carried out a variational comparison 
between the two superposition states 


I planar) 


1 

71 


(IQ®) + *IQy)) ; 


and 


|diag) 


1 

71 


(|Q,)+e*2V3|Q^)+g-i2V3|Q^)^. 


It is found that the latter has lower energy under Bogoliubov approximation. 
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One interesting consequence is that the angular momentum (L) order in the 3D p-orbital condensate state is 
noncollinear, which is different from the 2D case. The polarization configuration of (L) is shown in Fig. Such 
configuration exhibits non-zero chirality defined to be 

Xijk — Li • [Lj X -L/j), (45) 

where ijk denote nearby three sites of the four corners of a square plaquette in a clockwise direction. With ther¬ 


mal/quantum fiuctuations, the presence of such chiral order may lead to unconventional phase transitions (Li et al 
2014aD. 


In a relative shallow lattice, Eq. (42) derived under the harmonic approximation would receive significant unhar¬ 


monic corrections. There Gutzwilier calculations suggest a more exotic condensate with nematic order (Collin et al. 


2010), which spontaneously breaks the cubic lattice symmetry. 



FIG. 5 Noncollinear angular momentum order in a 3D p-orbital Bose-Einstein condensate in a cubic lattice (Cai et al. 2012b) 
Red arrows indicate the polarization direction of the angular momentum (L). 


5. Renormalization group analysis 

Eluctuation effects on p-band condensates beyond mean field theories are studied with perturbative one-loop anal¬ 


ysis (Liu et al. 2013), where the partition function takes the form 


with 


+ [ W<C(4)C(3)</>a(2)^a(l) 

+ 52</>;(4)<^:(3)(A.(2)(/.^(1) 

+ 53[<^:(4)(/.:(3)<^,(2)<^,(1) + /i.c.]}. 


(46) 


(47) 


/ ^:^^(27r)^(5(k4 + k3-k2-ki)(5(w4 + W 3 -W 2 -Wi) and </>a(j) denotes </>„(wj, kj). For general 
lattices lacking of C 4 rotational symmetry, the energy potential parameters are not equal, Vx ^ Performing a 
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momentum-shell renormalization group (RG) analysis, the fields are split into fast and slow parts, 
and k)||k|<A/s. Following the standard Wilsonian RG procedure (integrating out fast modes and rescaling the 

effective action for the slow modes), the RG flow equations (or /3-functions) for the potential parameters to one-loop 

order are obtained to be 


uT 

= 2f^ + 45i©(A - 1/2) + h^{fy - 1/2), 

(IT 

= 2fy + Agfeify - 1/2) + g2Q{fx - 1/2). (48) 

Here the dimensionless parameters are defined as f = rm/K^ and g = gm/{27:) and 0(x) is the Heavyside step 
function. The RG flow equations for the quartic couplings are 

25f - 2^1, 

25f-2ff|, 

gl, 

253(5?+52)- ( 49 ) 

With bare repulsive interaction, these quartic couplings are all marginally irrelevant. However they could strongly 
modify the RG flow of before they renormalize to zero. 

In the region with fx(0) > | and fy(0) > | the solutions are 


dgf 

dl 

dl 

dh 

dl 

dh 

dl 


rx{l) 


= e 


2 i 


A(0 


f^a:(0) + f dl'e ^q4g?(l') + 32 (^ 0 ) 
Jo 

^»/(0) + [ dl'e ^*( 4 ^ 2 ( 1 ')+ 52 (^ 0 ) 
Jo 


(50) 


In this region, and fy quickly run to positive infinity. In the region with A < 5 and fy < ^, the solutions are 

fxil) = rx(0)e'^\ 

fy{l) = fy{0)e^\ (51) 


from which the behaviors of RG flow are also fully determined by initial values of In other regions, one-loop 

corrections play more important roles in making the eventual values of positive or negative. Numerical studies 
have found interesting regions in the phase digram where fa,(0) < 0 and f^(0) > 0 (or vice versa) flow to +oo 

and Vy -A -hoc. Depending on the flow directions of four states can be identified: (1) Complex BEG {r^ +oc, 

Vy -A -hoc) ; (2) px BEG {vx +oc, Vy -A —oc); (3) Py BEG ( Vx —oc, Vy -A +oc); and (4) vacuum (f^, —oc, 
ry -A —oc). 

The RG study sketched above does not really capture the TSOC state because g^ flows to 0, making an illusion that 
quantum fluctuations wash away the phase locking between px and Py components. However this is not physically 
correct. A more careful RG study requires introducing U{1) and Z 2 order parameters to characterize the fluctuation 
effects in the TSOC state. 


B. Mott states, orbital exchange and frustration of bosons 


In the strongly interacting regime, bosons localize and form Mott insulator phases. Unlike the “featureless” 5 - 
band Mott insulators, the p-band Mott insulators have orbital degrees of freedom. Details of preparation of p-band 
Mott states including relaxation dynamics are studied in (Challis et al 2009). The orbital ordering is governed 


by the orbital exchange interactions which result from virtual boson tunnelings. Here we will derive the orbital 
super-exchange interactions and discuss the orbital frustrations on certain lattice geometries. 
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1. Mott states with filling factor larger than 1. 


The procedure to derive super-exchange interactions is to take the local terms as the leading part and the hopping 
terms as perturbation. Consider the two dimensional p-band Bose gas for example. The local interaction is given by 


(52) 


It can be verified that the angular momentum operator Lz commutes with the local interaction, i.e., 

Thus the eigenstates of the local interaction can be chosen as states with definite angular momentum. For filling 

2 

factor z/ > 1, the degenerate eigenstates with lowest energy = are 

|+) = ^|vac), 


|-) = 


m: 

y/i^. 


vac) 


where The states |+) and |—) have angular momentum and —respectively. On a square lattice, 

the tunneling Hamiltonian in the transformed basis reads 


Ft = ^ \T,AS^)bl 


h / 


h.c. 


-\-x 


with the matrices 


T{x) = 
T{y) = 


t\\—t± 

1 +Z L 


t\\-\-t± tl 


5 

2 

2 



t\\ -\-t± 

2 

t\\ -\-t± 


J- 


(53) 

(54) 

(55) 


2 2 

The low energy sub-space is spanned by the product states 

|{s(r)})^0rk(r)), 

where s{t) = ± and r runs over all lattice sites. All the states in this subspace have the same energy to leading 
order in U and there is thus a macroscopically huge degeneracy. The corrections due to the hopping term Ht lift the 
degeneracy. The first order corrections vanish because Ht does not connect any states in the low energy sub-space. 
The second order correction is calculated by the standard perturbation theory, 

|(m|F*|{s(r)})|2 


AF(|{5(r)})) = 


E- 


F(0)(|{s(r)}))-F(0)(|m))’ 


(56) 


where \m) is a higher energy state orthogonal to the product states |{5(r}), and is the leading order energy. 
Keeping only tunneling between nearest neighbors as in Eq. ([^, AE(|{s(r)})) simplifies to 


AF(|{5(r)})) = ^ AF(|5(r)s(r'))), 


(57) 


where r and r' are adjacent sites. Calculating the energy correction on a two-site state AE{\s{t)s{t'))) is straightfor¬ 
ward. The energy corrections are 

AF(| + +)) = AF(|--)) 

_ 3 (v\t\i+tj_\ 


4 ( -U 

AF(| + -)) = AF(|-+)) 


+ 


-uA + 1) j ’ 


3 f i^(i^ + i)|f|| + - t_i_\ 


-u 


+ 


-uA + 1) j ■ 


( 58 ) 
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Then the correction A£i(|{s(r)})) is given as 

A^(|{s(r)})) = J.s(r)s(r'), (59) 

(r,rO 


with 


J. 


+ 2 ) 

2(z/ + l) iT 


Including this correction into the Hamiltonian, we get 

AH = 

(r,r'> 


(60) 


where ay is defined to be a^ 


— 7/-1 


PLzP, with P a projection operator P = |+)(+| + |~)(~|- The orbital super¬ 
exchange makes the staggered angular momentum ordering energetically favorable. 

It should be emphasized here that the energy corrections in Eq. (59) actually do not depend on the orientation 
of the link r — r' and that the effective Hamiltonian in Eq. (60) is independent of lattice geometries. Considering 
p-band Mott insulators on a triangle lattice, the effective orbital model is geometrically frustrated making both of 
ferromagnetic and antiferromagnetic correlations suppressed. 

The above analysis holds in the deep lattice regime. Eor a relatively shallow lattice, the degeneracy in local Hilbert 


space could be lifted up (Collin et ai 2010). Treating such effects as perturbations, based on well established results 


m 


transverse field Ising models ( Sachdevj " 2011[ ) the staggered angular momentum order in the Mott state is expected 
to be stable when the perturbations are reasonably weak as compared to the super-exchange. But we would like to 
emphasize that the competition of charge (atom number for neutral atoms) and spin orders in the shallow lattice 
regime may alter the above speculation and lead to potentially rich physics. 


2. Mott state with filling factor 1 

Eor Mott states with filling factor = 1, the convenient basis to calculate the super-exchange interaction is the 
Px^ Py basis, rather than the Px ± ipy basis. The generic form of interaction in Eq. (20) is used here. Like deriving 
super-exchange for filling u > 1, we need to calculate the second order corrections of nearest neighbor product states— 
|1, 0; 1, 0), |0,1; 0,1), |0,1; 1,0), and |1, 0; 0,1), where a notation 




(61) 


[hl{v)r^ [bl{r)r^ [bUr')r^ [bl{r')y 

-y^rUxlrnylm'^lrriyl 


I vac) 


is adopted to save writing. The zeroth order energy of these four states is Ui. The higher energy virtual states that 
fftwillcoupletoare^(|2,0;0,0) + |0,2;0,0)), ^ (|2,0; 0,0) - |0, 2; 0,0)), |1,1;0,0), ^(|0,0;2,0) + |0,0;0,2)), 
;^(|0,0;2,0) — |0,0;0,2)), and |0,0;1,1), with corresponding energies 2Ui + U 3 , 2Ui — U 3 , Ui + 2 U 2 , 2Ui + U 3 , 
2Ui — Us and Ui 2 U 2 - Eor the link with r' = r -h x, the second order energy corrections are given by 


A^(|1,0;1,0)) - -2t| ’ 

AE(|0,1;0,1)) =0, 


A^(|0,1;1,0))=A^(|1,0;0,1)) = -^. 


(62) 


(Note that the transverse tunneling is neglected here, for the reason that the longitudinal tunneling is enough to lift 
the degeneracy and is significantly stronger than the transverse one.) Mapping Px and Py orbitals to the pseudo-spin 
1/2 states, a =t and respectively, the effective Hamiltonian on this link reads 


77x = + x) + Mz Wz{r) + o-^(r + x )], 


(63) 
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with Ji = [{Ui + Us)-^ + {Ui - Us)-^ - {2U2)-^], M, = -f [{Ui + Ua)-^ + {Ui - U^)-^] . Similarly, the effec- 
tive Hamiltonian for the link r' — r = ^ is obtained as 


Hy = Jiaz{r)az{r + y)- [^^(r) + C7^(r + y)]. (64) 

Then the total effective Hamiltonian for filling z/ = 1 on a square lattice is 

AH = ^ Ji [az{v)az{v + :r) + az{v)az{v + y)]. (65) 

r 


With Ui = 3 U 2 = 3 ^/ 3 , the coupling Ji is positive and the ground state has an antiferromagnetic ordering with 
alternating p-orbitals (see Fig. [^. For a one-dimensional lattice, makes Px orbitals favorable due to the effective 
Zeeman term Mz (Eq. (lU)). We mention here that including the transverse tunneling would give rise to even richer 


physics, e.g., an XYZ quantum Heisenberg model can emerge ( Pinheiro et al] 2Q13[ ). 


One key difference between filling u = 1 and higher fillings is that the super-exchange interaction depends on the 
orientation of the link r' — r, which makes the orbital frustration on triangle/Kagome lattices even more interesting. 



FIG. 6 Illustration of the alternating px/py orbital order (Li et al. 2011b) for a p-band Mott insulator at unit filling. 


3. Phase diagram of p-band Bose-Hubbard model 


The phase diagram of p-band Bose-Hubbard model in a two dimensional square lattice is studied by quantum (Hebert 


et al. 2013) and classical (Li et a/., 2014a) Monte-Carlo simulations. For filling of two particles per site or higher. 


a second order quantum phase transition from antiferromagnetic Mott state to the TSOC state is found at zero 


temperature in the quantum Monte-Carlo study as well as in Gutzwiller approach (Collin et al. 2010 Larson et al. 


2009 Martikainen and Larson 2012). In the weakly interacting regime at finite temperature, the fluctuations are 


modeled by a phase-only model studied by the classical Monte-Carlo. It is found that the TSOC state develops a 
two-step phase transition to the normal state, a Kosterlitz-Thouless transition followed by a higher temperature Ising 
transition. Sandwiched between the two transitions is a time-reversal symmetry breaking non-superfluid intermediate 
state. By combining the numerical results from Monte Carlo in the weak coupling regime and the analytical exact 


result from mapping the Mott limit of the p-band model to the orbital equivalent of the Onsager Ising model (Onsager 
1944), the phase diagram for p-band Bose-Hubbard model in two dimensions is proposed (Fig. [^. We would like to 


mention here that strong correlation effects may give rise to exotic intermediate phases between p-band Mott insulator 


and superfluid states (Xu and Fisher 2007). 
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Chiral Mott 
Insulator 



FIG. 7 The schematic phase diagram of the two dimensional p-band Bose-Hubbard model with Filing factor < 2. The 
Chiral Mott and superfluid states have staggered angular momentum ordering. At zero temperature there is a quantum phase 
transition between the chiral Mott and superfluid states. At finite temperature, there is a chiral Bose liquid state which has 
angular momentum order but no superfluidity. Upon heating, the chiral superfluid undergoes a Kosterlitz-Thouless transition 
into the chiral Bose liquid, which subsequently undergoes an Ising transition at a higher temperature into a normal Bose 


liquid (Li et al. 2014a). 


C. Interacting p-orbital fermions 

1. Nested Fermi surface—FFLO state 


Searches for superconducting Fulde-Ferrell-Larkin-Ovchinnikov (FFLO) phases with spatially varying order pa¬ 
rameters have been attracting tremendous interest in both atomic gases and electronic materials. It appears that 
the parameter window for this novel state to occur in two or three dimensions is quite narrow for conventional set¬ 
tings (Radzihovsky and Sheehy 2Q1Q|). (In an optic al lattice, the spin imbalance window to reach the FFLO state 


larger and much progress has been made to find FFLO phases 


2007 


2007 


|Casula et al. 2008 

\ Feiguin and Heidrich-Meisner 

1 |2009t |Liao et al.\ \ 

2010| Orsol |2007t Parish et al. 


|Yang[ 2001). Nonetheless, long range order is prohibited due to fluctuations in one dimension, which makes it 


challenging to observe FFLO states even in one dimension. 

Due to intrinsic anisotropy of p-orbital wavefunctions, FFLO states in p-orbital fermion systems are found to occur 


in a wide window in two dimensions or even in three dimensions (Cai et al. 2011), not restricted to half Ailing. Here, 


we shall focus on a two dimensional square lattice, where the tunneling Hamiltonian is (Cai et al. 2011) 


■^0 ~ ^11 \_^\c,cr,Y^x,(j, 


r+ax Y C^,cr,r^?/,cr,r+ay + 


-//^n^(r) - - (n-f(r) -n;(r)) 


( 66 ) 


with Cx^a the fermionic annihilation operators for px (py) orbital with pseudo-spin a ^ and ria = a^a,a 

the density operators for spin a. The transverse tunneling t±_ (^ ty) is neglected for simplicity and in presence of 
spin-imbalance, this leads to perfect nesting of p-orbital Fermi surfaces (see Fig. [^. 

In atomic gases the pseudo-spin components are hyperflne states. The interactions between them maybe engineered 
by s-wave Feshbach Resonance. With the Feshbach Resonance, a matured technique in experiments, the induced 
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FIG. 8 Nesting of p-orbital Fermi surfaces on a square lattice (Cai et al. 2011 | . The Fermi surfaces of px (py) orbitals are 
vertical (horizontal) lines. Solid and dashed lines are for the minority and majority hyperhne states, respectively. Fermi 
surfaces are perfectly matched as 0, with pairing momenta ±Q = zt{Skf, Skf), with Skf = kff — kf^. By lattice 

rotation symmetry, Q' = {Skf, —Skf) is the other choice of pairing momentum. 


interactions are given as (Zhang et al. 2QlQa[ ) 

-^int ~ U ^ ^ \p^x ^T 

r 

^ ^ J ^x,Y ’ ^y,Y ~^^x^Y^y 

r ^ 

+ ^ ^xy + h.C. 


(67) 


At tree level, J, Vxy and U are related: J = ^, Vxy = With attractive interaction, [/ < 0, the induced Cooper 
parings are ^ — \c'x^^y^x\.^y') i ^y,Y — {^y^,Y^y^,Y}^ ^xy,Y — {^x^,Y^y^,Y}•> Q'lid ^yx^r — {^y^,Y^x].,Y}‘ The BCS mean field 
Hamiltonian is 


Hm = Ho 

Y 

r 

-^^yx,YCyf^r^xf^r + ^xy,Ycl.f.^Cyf.^ + /l.C.j 
- Vxy ^ + ^x,YCyf^r^li^r + * 

r 

The last term Vxy < 0 locks the phase difference between and Ay at 0, which plays a central role in making the 
superconducting coherence two-dimensional. Solving the gap equation numerically yields 

A^y,Ay^^0, (69) 

= |A|cos(Q-r). (70) 


The mean field phase diagram has been mapped out in Ref. ( Cai et ^ 2011). The parameter regime supporting 
FFLO states is considerably large in spin imbalanced p-band fermions. 
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2. Nested Fermi surface—density stripes 


Besides the superconducting stripes in FFLO states, p-orbital fermions also naturally support density stripe orders, 
again from Fermi surface nest ing (Fig. [^. Stripe and checkerboard orders are found in a system of spinless fermions 


loaded up to p-orbital bands (Zhang et al. 2Q12| ), where the Hamiltonian takes a simple form 


ra/3 

/i ^ ^ 9 ^ ^ 

ra r 


with tai 3 = [t\\^ap — Fermi surface nesting of this system is pictorially illustrated in Fig. Fermi 

surfaces of Px and Py bands are approximately perpendicular to each other, which greatly suppresses the Cooper 
instability. The reason is that for the spinless case, the onsite interaction can lead to only cross-orbital Cooper pairing 
that is antisymmetric in px and py orbitals. The nearly orthogonal geometry of the two Fermi surfaces now makes it 
impossible to condense such Cooper pairs at a single center of mass momentum. In contrast, each Px (py) particle-hole 
pair in the density channel composed of one particle and one hole within the Px (py) band benefits from the fermi 
surface nesting. To simultaneously condense particle-hole pairs in each orbital band, the wave-vector 


Qi 2 = i2A;F) 

is most favorable (see Fig. [^. 

To characterize the Fermi surface nesting effect observed in Fig.[^ one can look at the density-density correlations, 
which can be calculated by the field theory with partition function Z = Tre“^^ — f ^ ^)) exp(—S' f), 

and the action 

Sf = f dT'^'lp^{r,T){dr - IJ.)lpa{r,T) 

r,a 

+ y^W(C(r + a/3,r)t/)„(r,r) + /i.c.) 

ycx(3 

+5 'r)V’^(r, T)'ipy{T, T)V’a,(r, t). (72) 

r 

The density fields are defined as pQ,(r,r) = '0*(r,T)'0(r,r), and density-density correlations are given by 

T 

na/3(^) = {pa,qPl3,q) ^ 


with q = (q, iuj)^ 


Paq — 


Yi j drpir,- 




and ipa{k) = f dr t/’a(r, r)e jg uggfuj ^ split into two channels—number density (p+) and orbital 

density (p-): p±{<l) = px,q ± Py,q- The correlations in these two separate channels are defined as 


n± = j^{p±{q)p±{-q))- 


Summing up ring diagrams under random phase approximation (RPA), the correlations are obtained as (Zhang et al. 


2012 ) 


2y° 

n±(Q,0 = 

1 ± .a V 


(73) 


with given by 


xO = D(Ef)ln(^^), 
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in the limit of t±_ 0. Here D{Ef) is the density of states near the Fermi surface and ujd is some energy cutoff in 

the field theory. Due to the logarithmic divergence in any arbitrarily weak attractive (repulsive) interaction ^ < 0 
{g > 0) induce divergence of n+ (n_) at sufficiently low temperature. The divergence of n+ and n_ indicates long 
range ordering of charge density wave (CDW) and orbital density wave (ODW), respectively. Transitions to these 
density waves are studied with mean field theory, where the Hamiltonian is approximated by 

HmF = (c'a^v+ap + “ M X] 

ra/3 ra 

Eg ^ ^ {'^x,rEIy^r T '^y,rMx,r EIx^rMy,r) 5 

r 

with Ma,r = (^a,r)- Sclf-consistcnt mean field calculations confirm that repulsive and attractive interactions favor 
CDW and ODW, respectively. The density patterns of these density waves are shown in Fig. 



FIG. 9 Fermi surface nesting and density waves of spinless fermions on p-orbital bands (Zhang et al. 2012). (a) shows the 


Fermi surface nesting. Red (dark gray) and green (light gray) solid curves indicate Fermi surfaces of px and py orbital bands, 
respectively. The solid arrow shows the ( 2 / cf , 2 / cf ) momentum of particle-hole pairing simultaneously satisfying the nesting 
condition for both px and py bands, (b) shows the checkerboard density pattern at half filling, (c) shows the density pattern 
of the striped CDW/ODW phase lower than half-filling. 


The order parameter for the charge density wave phase is introduced by 

p(r) = + C.C.] + const, 

where <pi and 02 are complex valued fields slowly varying in space. The phenomenological free energy reads 


F = 

(fv (K\V(t) 

i=i,2 ^ 


/ 


■j\ +r\(t)j\ +u\(t)j\ 




( 75 ) 
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Here, incommensurate filling is assumed and the theory has an emergent U{1) x U{1) symmetry; otherwise the terms 
such as (0^ + c.c.) are allowed and the theory has lower symmetry. 

The effective couplings, r, and v in Eq. (75) have been connected to microscopic parameters by field theory 
calculations ( [Zhang et al. 2012). At low temperature, we have (r < 0, rt > 2 v)^ and the system is in a striped 


CDW phase with the wavevector or Q 2 spontaneously chosen. Assuming is spontaneously chosen there is an 
algebraic long range order in 0i, i.e., (0J(r)0i(r')) oc At higher temperature it is found that the striped CDW 

phase first melts to a nematic phase through an Ising transition and then to normal through a Kosterlitz-Thouless 
transition. 


3. Strongly correlated orbital models 


At half filling p-orbital fermions described by the Hamiltonian in Eq. ( 

71 

) exhibits a Mott transition with strong 

repulsive interaction, which is studied in Ref. (Wu 

to 

0 

0 

00 

cr 

Zhao and Liu 

2008). In the fermionic Mott state, like in 


the bosonic case, fermions are localized on each lattice site. As a result, the low energy physics is described by an 
effective model of super-exchange interactions. 

Considering a link <r, r’> the super-exchange interactions are determined by energy corrections on the states 
|1, 0; 1, 0), |0,1; 0,1), |0,1; 1, 0) and |1,0; 0,1), where a notation is taken from Eq. (61) with the bosonic operators ba 
replaced by fermionic ones Cq,. Suppose this link is in the x direction, the tunneling is then + h.c.^ 

with the transverse tunneling neglected. Erom standard perturbation theory, the energy corrections due to virtual 
fermion fluctuations are 


ae;(|0,1;0,i)) = ae;(|i,0;1,0)) = 0, 

^2 

Ai;(|l,0;0,l)) = AE(|0,1;1,0)) = 

Mapping px (py) to pseudo-spin t (i) states, the super-exchange interactions are given in a compact form as 

Ks = Jz(r^{r)a^{r') + const, 

with Jz = Rotating p-orbitals by an angle 6^, we have the following transformation 


with 


U{0) = 


^U{ 0 ) 


cos 0 sin 0 
— sin 61 cos 61 


(76) 

(77) 

(78) 

(79) 


Eor a link oriented at an angle 0 with respect the x axis, the super-exchange interaction reads as 

h-ts = Jz^z{r)^z{r'), 

with az = {0)azU{0) = sin(26l)cra, + cos{20)az. 

On a square lattice, the Hamiltonian describing the orbital order is 

Ks = Jz'^crz{r)crz{r + a-cc) + crz{r)crz{r + ay), 


(80) 


(81) 


which has the same form as the p-band Mott insulator of bosons. This Hamiltonian supports an alternating PxjVy 
order as shown in Eig. On a honeycomb lattice, the super-exchange Hamiltonian is 




(82) 


with 




rr - ^ 

7i — + -jCFz, 

^2 — + 2 ^^' 


73 = -(^z- 


( 83 ) 
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Here the summation iiicludes one set of the ‘A’ sublattices (Fig. 10). This model, dubbed quantum 120° 
model ( Zhao and Liu[ 2008 Wu 2008b and Wu et al 2007), is geometrically frustrated. The complication of 


this model originates precisely from the spatial nature of orbitals, which makes orbital degrees of freedom drastically 
different from real spins. 

In three dimensions on a diamond lattice, the p-orbital exchange interaction leads to an exact orbital Coulomb 
phase characterized by ice rules and emergent gauge structures ( jChern and Wu 


2011 ). 



FIG. 10 Illustration of quantum 120° model on a honeycomb lattice. Ti ,T 2 and T 3 denote three different super-exchange 
interactions. 


4. Anti-Ferromagnetic phases of spinor p-orbital fermions 


As motivated by understanding the role of magnetism in high temperature superconductors, studies of antiferro¬ 
magnetic transitions in s-band fermions attracted tremendous interest, but the transition temperature is still out of 
reach for current cooling techniques. One way to improve the transition temperature could be provided by considering 
p-band fermions. The antiferromagnetic transition of spin-1/2 fermions loaded in p-bands of a 3D cubic lattice are 
studied in (Wu and Zhai 2008), where half filling (three fermions per lattice site) is assumed. The Hamiltonian 
describing such a system is H = Hq + with 


Ho — ^ 




ra^ 


and 


H-.r, 


E 


u ^ ^ ^r,Ck:,t^r,Q:,4 T bF ^ ^ 




t 

r,a,t 


C 


t 

r,/34 


Cr, 





With strong repulsion, we can project to the low-energy subspace determined by i^int; projecting out high energy 
subspace will contribute to super-exchange interactions. The low-energy states are the four degenerate components 
of total spin-3/2: | ttt), (I itt) + I tit) + I ttt)) ^ 75 (I tit) + I tit) + I ttt)) , and | ttt), where a notation 

ased. This is manifestation of the Hund’s rule. It is quite involved to perform the 
quantum mechanical second order perturbation theory here. A more elegant way is to take the Brillouin-Wigner 
approximation where the super-exchange interaction is given by 

Hj = -PcHoPeHrJp^HoPG. (84) 


Here Pe and Pq mean projections onto excited and low-energy subspaces, respectively. 

The calculations are greatly simplified by two observations (Wu and Zhai 2008). Firstly, the hopping processes only 
take place within the same orbital. Secondly, all terms in Ho acting on the low-energy subspace create eigenstates of 
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i^int with the same excitation energy U + 2W. Then the super-exchange Hamiltonian Hj on a link < r, r' = r + x > 
is given by 


E 


U + 2W 







By symmetrizing cj q, ^ one can show 





Then the super-exchange Hamiltonian is obtained to be 

Hj = J (85) 

<r,r'> 

with Sr = \Pg Sass' cl^a,s^ss'Cr,a,s'PG, and the effective coupling 

J = 4(t| + 2ti)/(9C/ + 18VF) >0. (86) 

The effective description is an isotropic spin-3/2 Heisenberg model. We remind the reader that this is the model for 
half filling, with the full Hilbert space of each site being spanned by three p-orbitals and two spins. Hand’s rule reduces 
the low energy subspace to the total spin-3/2 space. The ground state of the system thus has an antiferromagnetic 
long range order. This antiferromagnetic order is destroyed by thermal fiuctuations when the temperature is above 
Neel temperature ^ J. 


D. Topological bands and nontrivial orbital states 


In optical lattice experiments considerable efforts have been made to create topological bands. Neutral atoms 
loaded in such bands would experience effective magnetic fields due to non-trivial Berry curvatures. These experi¬ 
mental developments are motivated by consideration of novel quantum many-body states such as quantum Hall states 
and topological insulators/superconductors. While previous experiments largely focused on manipulating different 
hyperfine states of atoms with synthetic gauge fields, recent theoretical studies (Dutta et al 2Q14a|b Li et al 2013 


Liu et al. 2Q16[ [2M4{ 2Q1Q| [^n et aL\ 2Q11[ [^12b Yin et al. 2015) point to alternate ways to achieve topological 


bands by considering high orbital states in the optical lattices of non-standard geometry. 


1. Topological sp-orbital ladder 


A one dimensional ladder composed of two chains of s and p-orbitals is shown in Fig. The two orbitals are 
level in energy, and other lower orbitals are energetically separated with a large gap, and thus can be neglected when 
considering the sp ladder. The Hamiltonian describing this orbital ladder system is given by 


Ho 




ts 


c.+i + h.c. - E,- 


(87) 


where Cj = [a^j),(j)] , with a|(j) and Op^ij) being fermion creation operators for the s- and pa;-orbitals on the 
A and B chain respectively. The relative sign of the hopping amplitudes is fixed by parity symmetry of the s and Px 
orbital wave functions. As depicted in Fig. EH the hopping pattern plays a central role in producing a topological 
phase. With a proper global gauge choice, tg, tp and tsp are all positive. Focusing on half filling with chemical potential 
/i = 0, the Hamiltonian is particle-hole symmetric under transformation Cj (—l)-^Cj. Heuristically, topologically 
non-trivial band structure of the sp-orbital ladder may be speculated by rewriting the staggered quantum tunneling 
as 


tsp EN {-i<Jy) Cj+i + h.c. 
3 






















26 


resembling the spin-orbit interactions when the s and p orbitals are mapped to pseudo-spin (1/2) states. The physics 
of the sp orbital ladder is also connected to the more familiar frustrated ladder with magnetic tt-Hux, but the sp-ladder 
appears much easier to realize in optical lattice experiments. 

In the momentum space, the Hamiltonian takes a suggestive form 

^{k) = ho{k)l h{k) • (j, ( 88 ) 

where ho{k) = {tp — tg) cos(/c), hx = 0, hy{k) = 2tspSm{k) and hz{k) = —{tp + tg) cos{k). Here, IL is the unit matrix, 
and ax, ay and az are Pauli matrices in the two-dimensional orbital space. The energy spectrum consists of two 
branches. 


E±{k) = hoik) ± ^hlik) + hlik). 


An interesting limit is that when tg = tp = tgp, the two bands are both completely flat. As the momentum k is varied 
from —TT through 0 to + 7r, crossing the entire Brillouin zone, the direction of the vector h{k) winds an angle of 27r. 


In the notation of Ref. (Wen, 2012), the sp-orbital ladder belongs to the symmetry group T, C), as it has 


both particle-hole and time-reversal symmetries in addition to the usual charge U{1) symmetry. At half filling, it is 
characterized by an integer topological invariant, in this case the winding number 1. 

A manifestation of the nontrivial band topology is existence of edge states. It is easiest to show the edge states 
in the fiat band limit, tg = tp = tgp = t, by introducing auxiliary operators, 0±(j) = [cip{j) ± ag{j)]/^/2. Then the 
Hamiltonian only contains coupling between and of nearest neighbors, 

Hq 2t 0l(jf)0+(jf + 1) + h.c. 


One sees immediately that the edge operators 0+(l) and (l)-{Ng) are isolated from the bulk, i.e., decoupled from 
the rest of the system. These modes describe two edge states at zero energy. Away from the fiat band limit, the 


wavefunctions of the edge states analytically constructed in Ref. (Li et ai 2013) are found not to confine strictly at 


the ends, but instead decay exponentially with a characteristic length scale 

C = 2/log {\i^/t^ + tgp)/i^/t^ - tsp)\) . 

Here, recall that the implicit length unit is the lattice constant along the ladder leg direction. For tgp = ^/tgtp^ which 
includes the fiat band limit, the decay length ^ vanishes and we have sharply confined edge states. 

A topological phase transition to a trivial insulator state can be driven by inducing a coupling b etween s and v 


orbitals, AH = Ay CjayCj^ which can be engineered by rotating the atoms locally on each site (Gemelke et al 


2010). For the coupling strength Ay greater than some critical value A^, Berry phase vanishes and the systeni 
becomes a trivial band insulator, and the zero energy edge states disappear. Such a phase transition can be detected 
by measuring the density correlation between two ends in experiments. 

Regarding practical experimental realizations, careful treatments of band structures and Wannier functions are 


required as the details of tight binding models could receive significant corrections beyond harmonic approximations 


(Eq. ([2l)) (jCanczarek et al 

|20I4). One controllable way to couple s and p orbitals is to use a one dimensional shaking 

lattice (jDutta et al 

. 2015b 

|Lacki and Zakrzewski 2013 Przysiezna et al 2015 Sowinski 2012 Strater and Eckardt 

20I5| |Zhang et al , 

20151 Zhang and Zhou 2014), which has recently been realized in experiments (Fort et al 2011 

Khamehchi et al z 

!0I5 fNiu et al 2015 Parker et al 2013 Weinberg et al 2015). The other way to systematically 


control the sp-orbital coupling is to consider a noncentrosymmetric lattice where the coupling can be turned on and 


off by manipulating inversion symmetries (Liu et al 2016). 


2. Topological serriinnetal from mixing p and d orbitals 


We now turn to two dimensions and study how degeneracy of higher orbital bands may give rise to topological 

phases (Liu et al 

2010 Sun et al 

20I2b|). Consider a double-well optical lattice of the configuration shown in 

Fig. 12 ( 

Sun et al 2012b 

). By the space group symmetry {D4) of the lattice, the two p-orbital states {px and Py) 


are degenerate with the lowest d-orbital (i.e., dx‘^-y^) at high symmetry points in the momentum space. The lattice 
configuration is found to exhibit degenerate p and d orbitals. Considering a square lattice with three orbitals on each 
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FIG. 11 
A chain 


2013). This ladder consists of two chains, A and B. The s orbitals of 


A one dimensional sp-orbital ladder ( Li et al] 
are level with the px orbitals of B chain. The inter-orbital tunneling has a staggering sign as shown. 


site {px^ Py and dx-^-y^)^ the Hamiltonian of the tight binding model takes the following form 

Ho=S'^ dldr - tdd T] (4^^r+a, + 4c^r+a„ + h.C.) 

r r 

+tl X! (d,rP®,r+a, + 4,rPy,r+a^ + h.C.) 

r 

^ ^ iplc^rPx.Y+aiy + P\j^YPy,Y+aix h.C.^ 

Y 

E {dlpx ,r+a^ pl;^Y^r-\-a.x 

Y 

+4P?/,r+a^ - Pl^Y^r+^y + ^'C*) , (89) 

where ax (a^) is the lattice vector in x {y) direction, and Px,y^ Py,Y and are fermionic annihilation operators for 
Px^ Py and dx‘ 2 -y 2 orbitals at site r. The amplitudes of tunneling between these orbitals at nearby sites are tddi t\\i 
t_L, and tpd- With a proper gauge choice, these tunneling amplitudes are all positive. Here a point group 1^4 and 
time-reversal symmetries have been assumed. This tight binding Hamiltonian can be realized by a double-well optical 
lattice potential 


(Sun et al 2012b 


V{x^y) = —Vi[cos{kx) + cos(%)] 

-^V 2 [cos{kx + ky) + cos{kx ■ 


ky)]- 


(90) 


A typical configuration and the experimental protocol to realize it are shown in Fig. 12 (Sun et al. 2012b). By the 


point group symmetry (D^) of the lattice, the two p-orbital states {px and Py) are degenerate at high symmetry points 
in the momentum space. By dialing the relative strength of Vi and V 2 , the two p-orbitals may be tuned to degeneracy 
with the lowest d-orbital (i.e., dx 2 -y 2 ). That corresponds to the control of the value of the band gap 6 in Eq. (89). 


Band structure calculation has confirmed that the relevant physics is captured by the tight binding model (Sun et al 


2012b). 


In the momentum space, the tight binding Hamiltonian becomes, 


0 ^ 

II 

M 

4,k>4,k)^(k) ^ 

du \ 

Px,]<. 1 5 

k 

Py,^ ) 

2tdd{cos kx + cos ky) 6 

2 itpd sin fca; 

2itp 

—2itpdsmkx 2t|| 

cos kx — 2t± cos ky 


-2itpd sin kx 

0 

2t|| cos ky 


( 91 ) 


i sin ky 
0 


with H(k) given by 
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Depending on the value of the energy difference S, there are two types of band structures for this model. For 
S > Atdd + 2 t|| — dx 2 _y 2 orbitals are weakly hybridized with p-orbitals; for 0 < (5 < Atdd + — 2tx, the orbitals 

are strongly hybridized. For the latter case, a band touching point between the top and middle bands shows up at 
{kx =0,ky = 0) (F point). This band touching point has non-trivial topological property, which is characterized by 
the Berry flux defined as the contour integral of the Berry connection in the momentum space. 


7n = ^ dk • An(k), 

with n the band index, C a close contour enclosing the band-touching point, and the Berry connection An(k) = 
i{u\^\d\^\u\^), where |rtk) is the eigenstate of the Hamiltonian H(k). The Berry flux 7 ^ is quantized to an integer 
multiplied by 27r, and only two cases 7 ^ = 0 or tt are distinguishable without any symmetry requirement due to the 
gauge choice in |rt^(k)). However, with space-inversion symmetry, we can restrict I\un(k)) = |i 4 n(~k)), with / the 
space-inversion operator. The Berry flux then becomes well defined up to mod Att ( Sun et al] 2Q12b| ). For the band 
touching point considered here, 7 ^ is 27r, and this band touching is topologically protected (in presence of symmetry). 
Filling fermions up to such a touching point gives rise to a topological semimetal. 

A more illuminating way to show the topological protection is to construct an effective two band Hamiltonian in 
the vicinity of F point. Near this point, the da, 2_^2 orbital band is far below in energy and can thus be eliminated. 


With standard perturbation theory, the effective Hamiltonian is given to second order as (Sun et al. 2012b) 


^ _ f 'H22 ^23 A ^ f ^21^12 ^21^13 

~ V ^32 H33 J - /i V ^31^12 nsiHis 


with fi the chemical potential of the topological semimetal. 
Hamiltonian takes the following form 


Further expanding momentum around 0, the effective 


'HeE = 


ti + ^2 




+ ^tskxkyax + 


K)t 
h -12 


{kl - kl)a^, 


(92) 


where ti = t\\ -\- 


u. 


pd 


—T, h = -tj_, and ts = 


2t: 


pd 


The absence of 


2 tn-2t^+4Ua-S - component is protected 

The energy gap near F point is 2|h|, with h a planar vector 


2t|| —2t_L-\-Atdd—5 

by time-reversal and space-inversion symmetries 
h(k) = {2tskxky, — ky)). The vector h forms a vortex configuration with winding number 2 in the momentum 


space. At the vortex core (the F point k = ( 0 , 0 )), it is guaranteed that h = 0 , which means the degeneracy (or band 
touching) point is topologically protected. 

A question naturally arising is whether the required time-reversal and space inversion symmetries can spontaneously 
break at low temperature. Renormalization group analysis (Sun et al 2Q12b[ |2QQ9| ) points to the spontaneous 
symmetry breaking of time-reversal, and a state with angular momentum order (ip^y + h.c.) is stabilized at low 
temperature, if the interaction is repulsive. Taking this into the effective Hamiltonian, a gap opens at F point. As a 
result, the topological semimetal gives way to an insulator state at low temperature. This insulator is topologically 
non-trivial with finite Chern number. If the bare interaction is attractive, the renormalization equation shows that it 
flows to the fixed point of zero (usually called a marginally irrelevant term). In other words, the topological semimetal 
phase is stable against any attractive interaction in the perturbative renormalization group sense. 


3. Nearly flatbands with nontrivial topology 


In the system described by Eq. (91) at low temperature, the developed angular momentum order generates an 
additional coupling between two p orbitals. 


AH = y] iApl^^py^k + h.c., 


which breaks time-reversal symmetry and thus allows the Chern number to be non-trivial. With the parameter choice 
^ = —4:tdd + 2 t|| + A — 2t||A/(4t|| + A) and t± = t||A/(4t|| + A), the energies of the top band at F and M points are 
2011). Varying A with tdd = tpd = ^|| = t fixed, they found that the ratio of the bandwidth/band 


equal (Sun et al. 
gap is minimized 


TJ 20) at A/t = 2.8 for the top band. The top and bottom bands carry opposite Chern numbers 




























29 



(b) 



2012b| . (a), the experimental setup to realize the 


FIG. 12 Optical lattice realization of the topological semimetal (Sun et al. 
lattice potential in Eq. (90) for V 2 IV 1 > 1/2. The linear polarization of the incident monochromatic light beam (solid blue 
line) encloses an angle a with respect to the direction normal to the drawing plane, (b), the optical potential for Vi = 2.2Er, 
and V 2 = 3AEr. The darker (lighter) regions represent areas where the potential is low (high). The dashed line marks one 
unit cell of the lattice. 


±1, and are thus topologically non-trivial while the middle band is topologically trivial. Such nearly Hatbands with 
nontrivial topology mimic the Landau levels of 2D electron gas in strong magnetic fields. The flatness is crucial 
to reach fractional topological states in lattice models. Further numerical investigations have shown that fractional 


quantum Hall states are supported when the Hatbands are Hlled at certain fractional Hlings (Sheng et al. 2011 Wang 


et al 2011). 


E. Numerical calculations on lifetime and stability 


Earlier discussion of the lifetime of p-orbital BEG based on Fermi’s golden rule calculation largely relies on single¬ 
particle picture, and may underestimate many-body eHects. Numerical studies based upon Gross-Pitaevskii approach 
(Eq. ([M|)) ( [Martikainen 2011 Xu et al. 2013) indeed Hnd interesting phenomena beyond the scope of Fermi’s golden 
rule treatment. To study the TSOC superHuid in continuous space, a variational condensate wavefunction is taken. 


V;o(x) = 


K 


(93) 


with K the reciprocal lattice vectors, Xk and Ik the variational parameters, and (Qy) the minima for the Px 
(py) band. This wavefunction is superposed of two Bloch functions, and it breaks lattice translation symmetry, as 
required to describe the TSOC superHuid. The key features of TSOC superHuid which are time-reversal symmetry 
breaking and staggered orbital current, as predicted based on tight binding models, are conHrmed in the numerical 


calculations for continuous space (Xu et al. 2013). 


The stabilities of the TSOC superHuid state are investigated within the time-dependent GP equation. 


idr1p{'K^T) = 




V’(x,r). 


(94) 


Rewriting '0(r, r) into condensate and Huctuation parts, 

7/(x, r) = 7/o(x, r) + ^q(x, T)e'^‘^ + 1 ;* (x, 


the time-dependent GP equation determines the dynamics of Huctuations ( |Xu et al 2013) 


/ Mq(x,r) \ ^ / Wq(x,T) \ 

^^"Gq(x,r)) -'"^^‘»Gq(x,r)T 


( 95 ) 
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with 




^(q) = - 


( 'C(q) sV’o 

V 5^0^ 'C(q )) ’ 

(^V + zq) 


( 96 ) 


+ V{'k) + 2g\tl)o\ 


Note that the vector q is the lattice momentum after doubling periods to make q a good quantum number, and that 
rtq and r’q are periodic—r4q(x + 2a.x) = rtq(x + 2a^) = 'i^q(x), 'r’q(x + 2aa,) = 'r’q(x + 2a^) = r’q(x). The eigenvalues 
of cFzICci determine the Bogoliubov spectra, which are studied for square and checkerboard lattices. The fluctuations 
would grow in time if the eigenvalues are imaginary, leading to dynamical instability. This instability is cross checked 
by simulating real time dynamics in the continuous space where the optical lattice is treated exactly by a periodic 


potential (Xu et al. 2013), beyond the standard tight-binding model approximation. 


For a square lattice, the TSOC superfluid state is found to be dynamically unstable unless the interaction strength 
is extremely weak. In presence of dynamical instability, the lifetime of the TSOC superfluid state in a simple square 
lattice could be tens of milliseconds, rendering that such a state is experimentally unreachable for the simple square 
lattice. This conclusion is fully consistent with the early experimental finding of a relatively fast decay of the p- 
orbital atoms in a qnasi-lD lattice system ([Mhller et a/. , 2007). In contrast, for the checkerboard lattice as used 


in experiments (Olschlager et al. 


2011 


Wirth et al. 


2011), when the lattice is not too shallow and the interaction 


is not too strong, the TSOC superfluid state is shown numerically to be dynamically stable. This is consistent 
with the long lifetime as observed in experiments. Similar improvement with superlattices is also found in one 


dimension (Martikainen 2011). When the interaction is stronger than some critical value, the TSOC superfluid is 


no longer dynamically stable even for the checkerboard lattice. Based on the dynamical stability, a phase diagram is 


predicted in Ref. (Xu et al. 2013), which is consistent with experimental observations. 


Another way to understand the dynamical instability is to look at the energy cost for fluctuations i4q, 'r’q, which 


takes the following form (Wu and Niu 2001), 


SE^ = 


j (Wq(x),q(x)) /Cq 


Uq(x) 

^'q(x) 


(97) 


The fact that the eigenvalues of cTzlCa are imaginary implies the matrix /Cq is not positive definite (although the 
reverse may not be true), which means that the variational ansatz in Eq. (93) is not a stable saddle point of the GP 


energy functional. This in principle indicates tendency of forming some crystalline ordering ( Li et al\ 20IIa). 

The other type of instability is Landau instability for the reason that there are always Bogoliubov modes causing 
the free energy to be negative for p-orbital BEC, which means the state is a local saddle point that can decay into 
the lowest s-band. However this instability is less important than the dynamical instability within the lifetime of 
experiments. The time scale for Landau instability to destroy the p-orbital BEC is estimated to be 500ms while it 
is found to be around lOms in numerical simulations for dynamical instability. Although the p-orbital BEC is not 
strictly a metastable state due to Landau instability, it is fairly stable within the experimentally relevant time-scale. 
In the checkerboard lattice experiment (Wirth et al. 2011) where each lattice site actually represents an elongated 


tube in the third direction, the dynamical phenomena are even richer. Eor example, a collision process with two 
atoms decaying into the lowest band is allowed as the energy could be released to the kinetic motion in the third 


direction (Paul and Tiesinga 2013). 


The dynamical instability of excited band condensate in a double-well lattice has also been studied in detail, and 
the loop structure in Bogoliubov spectra is found to be correlated with the dynamical instability (Hui et al. 2012). 


IV. EXPERIMENTAL PROBES AND NOVEL LATTICES 

The theoretical discovery of richness of many-body physics with p-orbital atoms has motivated considerable ex¬ 
perimental efforts in recent years. So far the experiments have been done only for bosonic atoms. It has been 
demonstrated in a checkerboard optical lattice that the chiral p-\-ip Bose-Einstein condensate gives rise to nontrivial 
quantum interference. In this section, we will review the experimental challenges to detect the chiral order, the recent 
proposals in theory and attempts in experiment, and the current status. 
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A. Early experimental observations of higher bands in a cubic lattice 


Coherent bosonic cold atoms were observed in the higher bands of an optical lattice in the pioneering experiments 


of accelerating lattices (Browaeys et al. 2005) and of cross-band Raman transitions (Muller et al. 2007). 


In the experiment of Muller et al 2007 the sample is prepared by first loading a Bose-Einstein condensate of ^' Rb 


atoms into a deep symmetrically simple cubic 3D optical lattice formed by three far detuned laser standing waves. 
For this deep lattice, it can be treated as an array of 3D harmonic oscillators with discrete vibrational levels, which 
can be labeled as \mxmymz) with rrij the vibrational quantum number along the j axis. Population transfer in these 
orbital levels can be controlled using a stimulated two-photon Raman process with propagating laser beams along the 


X axis (see Fig. 13), which provides an inter-orbital coupling 




with fleff the effective Rabi frequency. The experiment restricts the Raman coupling to the lowest Bloch bands and 
demonstrates orbital transition from the |000) state (s-orbital) to |100) (pa,-orbital). Rabi oscillations between the 
two orbitals have been observed. A maximal transfer efficiency of nearly 80% is achieved. 

The decay of atoms into the lowest orbital due to collisional events has also been measured. The lifetime was 
found to be 10 — 100 times longer than the tunneling scale. Emergence of coherence compatible with a Bose-Einstein 
condensation to a nonzero momentum state has been seen; yet the experimental system was anisotropic and the 
predicted Px + ipy-wdive condensate was not studied for the absence of Px and Py orbital symmetry. 



(b) |ooo> |100> 




X 



^Raman pulse 


FIG. 13 Population of higher orbitals with Raman transition (Muller et al.\ 2007). (a), schematic of stimulated Raman 


transitions from s- to p-wave orbital, (b), the population of the lowest (i) and first excited band (ii) measured by time-of-flight 
techniques. Rabi oscillations between the s- and p-wave orbital demonstrate the coherent coupling. 


B. Observation of high-band condensation in a checkerboard lattice 


After the early observation of higher band population ( Johnson et al] 2QQ9[ [Muller et al] 2QQ7|), long-lived Bose- 


Einstein condensate in the high-bands was not achieved until the groundbreaking experiment ( jWirth et al. 2011). 


In this experiment, a square optical lattice, composed of two classes {A and B) of (tube-shaped) lattice sites is used 
(see Fig. 14). Formed by two standing waves oriented along the x and y axes with polarization along the ^ axis, the 
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lattice potential is 


V{x,y) = —^ \r][{zcos{a) + y + eze 

z{e^'^y + ee-^'^y)f , 


(98) 


where y ~ 0.95 accounts for a small difference in the powers directed to interferometer branches, e ~ 0.81 accounts 
for the imperfect retro-reflections, and the angle a permits tunability of anisotropy in the x-y plane. An isotropic 
p-band with degenerate band minima arises when cos(a) e (or a = o^iso ~ ^r/S). The controllability of the phase 
difference 0 allows to adjust the relative depth of potentials at A and B sites, which is crucial in this experiment to 
populate higher bands. For 0 < ir{2 the A sites are shallower than the B sites and vice versa. 

Initially a Bose-Einstein condensate of rubidium (^^Rb) atoms is prepared and the lattice potential is adiabatically 
turned on with 0 = 0.387r such that 5-sites are much lower than A. A lowest band lattice Bose-Einstein condensate is 
thus created with most of atoms confined in B sites. Then 0 is rapidly increased to a final value Of > 7r/2 such that the 
5-orbitals in the B sites are level with the p-band of the lattice in energy. In doing so, atoms are efficiently transfered 
to the p-band. Since this preparation procedure is abrupt, the prepared state is not immediately a condensate state 
but rather an incoherent state, in which the atomic distribution in the Brillouin zone is fairly uniform. Surprisingly, 
after some holding time around lOms, sharp peaks arise at p-band minima and the p-band Bose-Einstein condensate 
spontaneously emerges. In theory the emergence of phase coherence is beyond the scope of Gross-Pitaevskii approach, 
and can be studied by constructing a quantum rotor model (San et al 2012), where the dynamics is well captured 


by the truncated Wigner approximation (Polkovnikov et al. 2002). 


The p-band condensate is not a ground state of the system but a metastable state; decaying into the lowest band 
is unavoidable. In this checkerboard lattice, the band gap between p-band and the lowest band is largely mismatched 


with the gap between p-band and the higher band, and Eermi’s golden rule calculation (see Sec. II.C) predicts a 
significant improvement of stability. In the experiment, the lifetime of the p-band condensate could reach lOOms or 
longer. 

Erom the measurements of momentum distribution, the experimental evidence of p-band condensate is conclusive. 
However there is no direct evidence for the orbital ordering in the TSOC state as predicted in theory. As a step 
further, a phase diagram is mapped out with va rying a (controlling the anisotropy) and the phase diagram is quanti¬ 
tatively consistent with theoretical predictions (Olschlager et al. 2013). The remarkable consistency of experimental 
observations with theories strongly suggests the p-band condensate be a TSOC state. Yet, direct evidence of the 
orbital order requires further experimental investigation. 


Population of even higher bands, say /-bands, is also achieved in this checkerboard lattice (Olschlager et al. 2011) 


thanks to the tunability of relative depth between two sublattices. Similar procedure was implemented as in preparing 
the p-band condensate. The resulting /-band condensate also has a complex nature. The condensate wavefunction 
locally resembles the superposition ip[ 3 ^o] ± ^'0[o,3] of eigenfunctions '0[n,m] of a 2D harmonic oscillator with n and m 
oscillator quanta in x and y directions, which has a spatial {2x^ — 3x) ± i(2p^ — 3y) dependence locally. The complex 
/-band condensate emerges from the same mechanism as the TSOC state of the p-band, namely maximizing the local 
angular momentum. 

Besides the way of loading atoms into the excited bands demonstrated in the checkerboard lattice, there are other 


vibrating lattices (Lacki and Zakrzewski 2013 Sowinski 2012). 


possibilities, for example by Bloch oscillation techniques (Larson and Martikainen 2011 Tarruell et al. 2012) or by 


C. Early experimental realization of double-well lattices 


Observations of higher bands in optical lattices are achieved in the early experiments manipulating double-well 
lattices, which were largely motivated by implementing coherent control of quantum degrees of freedom (Anderlini 


et al. 2007 

Cheinet et al. 

to 

o 

o 

00 

1 Lundblad et al. 2008 

Sebby-Strabley et al. 2006 Trotzky et al. 2008) 

Here we use the experiment (I 

Sebby-Strabley et al. ' 

2006) to demonstrate how the higher bands are p 


double-well lattices and what consequent observables are achieved. This double-well lattice is a two dimensional lattice 
formed by superimposing two lattices with orthogonal polarizations. Having a laser setup as shown in Eig.jl^a), the 
electric field generated by the four laser beams is Re[5(x, with 


E{x, y)=E (e*''® + ei 


^ikx _j_ ^i{29-\-2(f)—kx) 

+5 


62, 


(99) 
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FIG. 14 Population of excited bands (figure provided by A. Hemmerich as courtesy), a, the checkerboard lattice with two 
sublattices A and B. b, the experimental sequence to populate excited bands versus the final value of 0 (see Eq. (98)) in step 
2 in b. c, the populations of higher bands with varying Of. The upper panel illustrates momentum distributions in different 
Brillouin zones (top row), and their dependence on d shows the condensation in the A+ point and the band relaxation to 
the hrst BZ after long times. In (e) three momentum spectra are shown, with the middle one c orresponding to the interesting 
case of equal populations in A+ and X-. Original results in a different form were published in (Wirth et al. 2011). 


where k = 27: jX (A is the wavelength of the laser light), 0 = kdi + 50^ and (p = kd 2 + Sp (the extra phase shifts SO and 
Sp are polarization dependent and can be controlled in experiments). We have neglected several imperfections such 
as imperfect alignment and reflections for simplicity here. In experiments these imperfections could cause technical 
challenges. For light polarizations being all in plane such that ei = ^, 62 = f, we have a laser intensity field 


y) /lxy,0 ( 100 ) 

= 2 cos{2kx — 20xy — 2(pxy) + 2 cos(2% + ‘2<pxy) + 4, 


with subscripts in 0 and p specifying the polarization dependence. For the out-of-plane case, ei = 62 = the laser 
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intensity field is 


h{x,y)lh^o ( 101 ) 



.k , , Oz 

2 

k Oz 

16 

cos{-{x + y) - — 


cos{-{x-y) - — - (l)z 


The laser field creates an optical potential V{x,y) oc {Ixy{x^y) + Iz{x^y)). With in-plane and out-plane polarized 
laser beams combined, a double-well lattice can be created (Fig. p!^b)). 

Ground state can be achieved by adiabatically loading atoms into the lattice. For the double-well lattice, different 
from simple Bravais lattices, the band gap could be very small compared with the energy scale with Tioad Ihe 

loading time. Then the Landau-Zenner transitions across the lowest and first excited bands can be significant. The 
population of the first excited band causes the oscillations in the momentum distribution measured in time-of-fiight, 
which are observed in experiments. 



FIG. 15 Laser beams to generate a double-well lattice ( Sebby-Strabley et al\ 2006| ). (a) shows the laser setup. The incoming 
beam with wave vector ki is reflected by mirrors Ml and M2 and after traveling distance di returns to the cloud with a wave 
vector k 2 . The beam is then retro-reflected by M3 and returns with a wave vector ka, having traveled with an additional 
distance 2^2. (b) shows the generated double-well lattice with G,o//xy,o = 0.4, 4>xy — (t>z — 7r/2 and Oxy —Oz — —7r/2 (see text). 
The darker (lighter) regions represent areas where the potential is low (high). 


The relation between the observed oscillations in the momentum distribution and the population of the excited 
band can be quantified by constructing a two-band model, 

H =J24>lTrr'4>r', (102) 

r r' 

with (/)r = [0A,r5 whcrc 4>a and <pB are annihilation operators for the localized orbitals, reA(x —r) and reB(x —r), 

in the two sub-wells at site r in the double-well lattice. In momentum space, the Hamiltonian then reads H = 
0'*'(k)H(k)(/)(k), with (/)(k) Fourier transform of 0r- After loading bosonic atoms into the lattice, the condensate 
is a superposition of the ground state and excited state at lattice momentum k = 0, 


U(0) = holL -h hxCFx + hydy, 


Writing H(0) as 
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the dynamics of the state j^/;) is given as + ipee with A = 2^h‘l + In terms of 

0 A,B(k) basis, we have 


IV’(i)) = { [^eS cos( 7/2) - sin(7/2) (/)Ji(0) 

+ sin(7/2) + 003(7/2) (/^t(0)} |0)> 

with 7 the polar angle of the vector (hz^hx). The momentum distribution is then given as 

n(k) = const + 2Re cos( 6 >/ 2 ) + '?h^(k) sin( 6 >/ 2 )) (reB(k) cos( 6 >/ 2 ) + re^(k) sin( 6 >/ 2 ))] 

where reA,B(k) is the Fourier transform of — r). The population fraction of the excited band could thus 

be extracted from the dynamical evolution of momentum distribution. 

Although the above discussions were restricted to the setup in the experiment (Sebby-Strabley et al, 2006), the 
coherent oscillation in time-of-flight is a generic phenomenon when a superposed state of ground and excited bands 


is prepared. And indeed similar oscillations are observed in other double-well lattices as well (Anderlini et al. 2007 


Muller et al, 2007 Trotzky et al 2008). 


D. Theoretical understanding of experiments 


Early theoretical studies of p-band condensates focus on the case with the point group D/^ symmetry. For the lattice 
potential realized in the experiment of Hamburg (Eq. ([9^), the point group symmetry is maintained only for the 


ideal case e = 1 and a = 0, where the potential reduces to E = — Vq {r]‘^ cos^ kx + cos^ ky ^ 2 r] cos 0 cos kx cos ky). 
For the realistic situation with e < 1, the D 4 symmetry is thus broken and only reflection symmetry with respect to 
the x-axis is preserved. The asymmetry could be partially compensated by setting a = aiso, for which the potential 
reads V = —Voe [ 7^6 cos^ kx + cos^ k y] — Vpeycoskx \cos{ ky + 6 >) + cos{ky — 0)] . The consequences of asymmetry 
are studied in detail in (Cai and Wu 2011 Shchesnovich[ 2Q12| ). 

The band structure is calculated by plane-wave expansion ( |Cai and Wn[ 2Q11| ). The reciprocal lattice lattice 
vectors are defined as = mbi -|-nb 2 , with bi ^2 = (±7r/a,7r/a) (a the lattice constant). Taking the single¬ 

particle Hamiltonian Hq = —k?V‘^/(2M) + E(x), the diagonal matrix elements are (k -1- G^^n|^o|k + G^^^) = 
Er {[akx/i^ + (m — n)]^ + [akyli: + (m + ^)]^} , with the single-photon recoil energy, and the off-diagonal matrix 
elements are 


(k|iJo|k + G±i,o) = -^7?e(cosae'F*® + 
(k|iJo|k + Go.±i) = -^7?(cos(Q;)e^*® + 

{k|iJo|k + G±l,:pl) = 

(k|ifo|k + G±i,±i) = -^ecosa. 


There are four time-reversal invariant points in the Brillouin zone, O = (0,0), X± = (±^, ^), and M = (^, ^), 
at which the Bloch functions are real valued. The band spectra are symmetric at these points, and consequently 
3k^(k) = 0, which means that they are saddle points in the band structure. For the choice a = aiso, the second 
band has double degenerate minima at and X_. For a < ai^o (<a > aiso) 5 {^-) becomes the unique band 

minimum. 

To investigate the interaction effects, the Gross-Pitaevskii equation 


' 2M 


Eeff(x)U(x) = F;T(x), 


(104) 


with Eeff(x) = E(x)+pp|T(x)p, is solved self-consistently by assuming the condensate wavefunction is a superposition 
of Bloch functions at X ±, 


T(x) = cos((5)'0x+(x) + sin((5)e*'^'0x_ (x). 


(105) 
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The Bloch functions 'iIjx± have nodal lines in space, while the variational condensate wavefunction could avoid nodal 
lines by having complex values (with (5 7 ^ 0 or 7 r/ 2 , and 0 7 ^ 0). The complex solution is spatially more uniform and 
thus more favorable by interactions, but at the same time costs more kinetic energy when a 7 ^ aiso- 

The competition between interactions and anisotropy leads to an interesting phase diagram containing two real and 
one complex states of Bose-Einstein condensation. The Gross-Pitaevskii approach finds second order transitions at 
zero temperature (Cai and Wu 2011). The phase transitions can be understood within a Ginzburg-Landau theory. 


F = -ri\^p+f - r2\^p-\‘^ + 9i\^p+\‘^ + g2\ip-\‘^ 

+93\^+f\'>P-f+94i^f^- +C.C.), (106) 


with describes the condensate component at X±. The Umklapp term ^4 > 0 favors the complex state. Assuming 
ri,r 2 , and — 2^4 > 0 , the complex state occurs in the regime 


gs - 2^4 ^ ^ 2^1 

2^2 r2 93 - ' 


The predicted phase diagram is confirmed in the experiment (Olschlager et al. 2013). 


(107) 


E. Measurement of orbital orders by quench dynamics 


Direct measurement of orbital ordering, namely the staggered angular momentum, was thought to be an experi¬ 
mental challenge, which motivates a theoretical proposal of using quench dynamics ^lT et aL \ |2014a| ) . The key idea 
could be understood by drawing an analogy between the two orbital states at each site {Px^Py)^ and a pseudospin- 1/2 
degrees of freedom (t, i). In this analogy, the Px ± ipy state corresponds to a pseudospin pointing along the y direc¬ 
tion in spin space. Applying a ‘magnetic field’ along the x direction to this pseudospin should then induce Larmor 
precession, leading to periodic oscillations of the ^-magnetization, corresponding to the population imbalance between 
two p-orbitals, AN = N{px) — N{py). Here we consider a square lattice. We can take a certain initial state and then 
quickly turn on a strong ‘magnetic field’ 


i^mag = y](-l)’'‘'+'’‘'A(r) [bl{r)by{r) + h.c.] 
r 


(108) 


at time r = 0. For simplicity, the ‘magnetic field’ is assumed to be strong enough to completely dominate the short- 
time dynamics. If initially a staggered superposition px ± e'^^py is prepared, all local Larmor precessions add up to 
produce a macroscopic oscillation in the orbital imbalance AN. This imbalance evolves within a Heisenberg picture 
as 


dAN{r, t) 
dr 


-i[AAf(r,T),iJi„ag] 

-2A(r)Lrg(r,r), 


(109) 


with the staggered angular momentum operator, whose time evolution is described by 

jrstag 

=2A(r)AAr(r,T). (110) 

dr 

This leads to oscillations in (AA^(r,r)), 


{AN{r,T)) 

= (AA7(r,0))cos(2A(r)T) - (Lf'‘S(r, 0)) sin(2A(r)T) 

= A{r) cos( 2 A(r)r + ( 111 ) 

where (AA(r, 0 )) and (I/|^^s(r, 0 )) denote the orbital imbalance and staggered angular momentum for the initial state. 
The trigonometric form of this time-dependent equation thus defines the quantities A(r) and 0(r), ready to compare 
with the experimental measurement of AN. 

Neglecting spatial inhomogeneity in A(r) and 0(r), we can set A(r) = A and (/)(r) = 0, and extract the initial 
angular momentum order from the amplitude A and the phase shift 0 in the dynamics of the spatially averaged orbital 
imbalance {AN{r)) = l/Ng AA(r, r). The coefficient A can be read off from the oscillation period tq = tt/A. The 
orbital population imbalance can be measured directly in time-of-flight experiments. 
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For a C 4 symmetric initial state with non-zero staggered angular momentum, but no orbital imbalance, {AN{t)) is 
expected to oscillate with a non-zero amplitude and phase shift 0 = ±7r/2 whose sign will fluctuate from realization 
to realization. By contrast {AN{r)) = 0 should be observed for a completely disordered state. The amplitude of the 
signal is thus a direct measure of the staggered angular momentum order parameter. 

In the case that C 4 symmetry is explicitly broken as achieved in the recent experiment, a state with an initial orbital 
imbalance but no angular momentum order would exhibit oscillations with a finite amplitude but no phase shift, i.e., 
0 = 0. In contrast, for a state with angular momentum order, the spontaneous time-reversal symmetry breaking 
yields a finite phase shift 0 7 ^ 0 , which would vanish in a singular fashion as we tune from the angular momentum 
ordered to disordered regime through a second order phase transition. 

The required coupling i^mag can be engineered by adding a quench potential Vkiag(x) modulated in the (1,1) 
direction with respect to the original lattice potential. The add-on potential generates a coupling between px and Py 
orbitals 


e(r) 


h 



3ix+3Ly 

) 

[ ] 

J 


dmcjo a‘^dP 


1 1 —^0 5 


( 112 ) 


where ujq is the harmonic oscillator frequency of the lattice wells hosting the p-orbitals and a = ja^,! = ja^j is the 
lattice constant. The above estimate for the coupling strength is valid in the tight binding regime when the quench 
potential is weak as compared with the original optical lattice. Without loss of generality, one may consider an add-on 
optical potential of the form 


14iag(x) = -T COS^ 


2 z/ + 1 


{Kx + K^) • X 


(113) 


with some integer z/ > 0 , a positive amplitude T, and (K^,, K^) denoting the primitive vectors of the reciprocal lattice. 
This potential leads to a PxjVy coupling 

= (114) 

with Ex the photon recoil energy with wave number l/2|Ka, + K^|. The staggering factor in the engineered coupling 
is crucial to probe the staggered angular momentum order. 

This quench proposal brings other interesting possibilities in addition to providing a method to probe orbital order. 
For instance, one can simulate spin dynamics in solid state materials by studying orbital dynamics of p-band bosons. 
One advantage about orbital dynamics is that engineering artificial effective magnetic fields is intrinsically easier due 
to the the spatial nature of orbital degrees of freedom than engineering real staggered magnetic fields. 


F. Measurement of the complex phase by Raman transitions 


There is another proposed scheme to measure the inter-orbital phase coherence in px ± ipy superfluid by Raman 
transition (|Cai et al. 2012a). In the Px Eipy superfluid, condensation takes place at and X_ and the condensate 

state is |T) ex f j |0) in general. The idea is to transform the phase coherence to number difference in 

momentum space. With a Raman operation, bosons in the original condensate can be transfered to a state with 


^'x+ = 
b'x- = 

With 0 = 0, the phase coherence in the Px ± ipy state is then transformed as 

{ib^xM. + h.c.) = {hPb'x_ - bU'xJ = Sn', 


(115) 


(116) 


which can be extracted in time-of-flight experiment. 

The required Raman transition can be implemented by two traveling-wave laser beams along different directions 
with corresponding wave vector ki ^2 and frequency uji ^2 ( Dnan[ 2006). These laser beams induce an effective Raman 
Rabi frequency with a spatially varying phase Q(x,t) = where = ki — k 2 , Scj = cui — CJ 2 , and 0 
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FIG. 16 Illustration of proposed Raman scheme to detect the complex orbital order in px di ipy superfluid ( |Cai et oT 2012a). 
(a) shows the Raman pulses with different propagating directions to build up momentum transfer between bosons at and 
X -. (b) and (c) show the time-of-flight imaging after Raman transition for the complex coherent px =b ipy state and incoherent 
mixing of px and py condensates, respectively. 


is the relative phase between the two laser beams (see Fig. [T^a)). The effective Hamiltonian for the Raman process 
is described by 


Hr = J (ixO(x, t)0^(x)0(x) + /i.c., (117) 

where (/)(x) is the boson annihilation operator in continuous space. The generated spatially dependent potential 
couples the two condensate components at the two momentum points [in the Hamburg experiment 
2Q11| ) X± = (±7r/2,7r/2), requiring — X_ = (tt, 0)]. 

To avoid complications of interband transitions (with band gap A) and dynamics caused by tunnelings (t), an 
optimal choice for the Raman coupling strength is t ^ HQo A. For the experimental situation, the Raman 
coupling strength should be chosen to be Qq ^ 27r x 0.5kHz. Thus the required duration of the Raman pulse is around 
1ms. To get efficient Raman operation, the frequency Suj should match the energy difference between the initial and 
final states which is around a few Hz. Therefore the phase accumulation Sut within the duration of Raman pulse is 
negligible. With this approximation the Raman coupling is simplified to be 


(Wirth et al 


+ h.c. (118) 

k 

Here A(k) is the k dependent effective coupling, which can be calculated from the Bloch functions. For the Hamburgh 
experiment, it is estimated that X{X±) ^ 0.98flo = A. Choosing the duration of the Raman pulse to be X6t = 7r/4, 
the required state transfer in Eq. ( |115[ ) is achieved. The resultant density difference is 

6n' = {ie'^^bj^_^bx_ + h.c.). (119) 

For the Px ± ipy superfluid, the density difference would be 6n' oc cos((/)). With 0 = 0, 6n' = {ib\^_^bx_ + h.c.) 
represents the order parameter of the complex orbital ordering (Fig, 
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G. Interference measurement of the complex phase 


In a recent experiment (Kock et al , 2015), that generalizes the idea of Young’s double slits, an interference measure¬ 
ment has been implemented to detect the inter-orbital phase coherence in the Px ^ipy superfluid. In this experiment, 
two independent copies of the lattice condensates are prepared with the experimental setup as illustrated in Fig. EZl 
The condensates are simultaneously prepared in the second band in two spatially separated regions of the lattice. 
After the state preparation, all potentials are switched off. The zeroth-order Bragg peaks observed in the x^-plane 
carry interference patterns in the 2 ; direction due to overlapping contributions from the condensates originally separate 
in space. In the simplified picture approximating the two condensates by two point sources, the wave length of the 
density grating in the interference is ’ with txoF the time of ballistic expansion, dz the spatial separation 

of the two condensates. This estimate is quantitatively consistent with experimental results. 

In the ballistic expansion, the Bragg peaks (labeled by 1, 2, 3 and 4 in Fig. [T7| yield the Fourier components of the 
condensate wavefunction, and we can associate a phase for each component, 6 >j=i^ 2 , 3 , 4 - Since the spatially separate 
condensates are decoupled, they carry different phases, Oj and O'y From the relative phase = Oj — Oy we can 
introduce AOij = AOi — AOj, which directly determines the correlation among the interference patterns in the Bragg 
peaks. If AOij = 0 (tt), the density patterns of the ith and jth peak are positively (negatively) correlated. The 
interference patterns obtained in experiments yield that AOi^s = ^^ 2,4 = 0, over 420 independent realizations, and 
that A0i^2 = A 6 >i ,4 = A 6 > 2,3 = ^^ 3 , 4 ^ their value spontaneously chooses 0 or tt. The interference measurement 
unambiguously tell that the phase of different momentum components is indeed correlated. To the best of our 


knowledge, the experiment (Kock et al. 2015) appears to be the first phase sensitive measurement which poses an 


important constraint on the nature of p-orbital Bose Einstein condensates. It is desirable that future experiments can 
directly probe the phase lock between the condensate components at two band minima, corresponding to the Xy and 


X- points in the paper (Kock et al. 2015). 


V. DISCUSSION AND OUTLOOK 


A. Orbital physics in electronic materials 


The crystal structure of the atomic ions in solids provide confining potential for electrons due to strong Coulomb 
force. Electrons in solids are usually nearly localized on atomic ions and the resulting orbital wavefunctions (or 
the shape of the electron cloud) are determined by the strong confining potential. This orbital degree of freedom 
is of great importance in correlated materials such as transition metal oxides ( Tokura and Nagaosa[ |2QQQ| ). Many 
intriguing phenomena such as metal-insulator transitions and colossal magnetoresistance can be attributed (or partially 
attributed) to the interplay of d-orbitals with charge and spin degrees of freedom. 

Considering a transition-metal oxide material with perovskite crystal structure, d-orbital electrons localized on the 
transition-metal atom are surrounded by six oxygen ions which give rise to crystal field and consequent energy 
splitting of the d-orbitals. Orbital wavefunctions pointing towards the negative-charged oxygen ions (the eg orbitals, 
dx 2 -y 2 and d^z^-r^) have higher energy compared with those pointing in other orientations (the t 2 d orbitals, dxy^ 
dyz and dxz) due to Coulomb repulsion (see Eig. 18). The spatial nature of orbital makes it intrinsically attached 


to the crystal fields, even in the absence of the relativistic spin-orbital interaction, and this intrinsic coupling of 
orbital degree of freedom to crystal fields and the resultant crystal symmetry make it distinct from real spins. When 
orbitals are modeled as pseudo-spins, the model Hamiltonian is in general lack of SU ( 2 ) symmetry. Consider a typical 
Mott insulator LaMnOs as an example. A neutral Mn atom has an electron configuration 3(i^4s^. Losing three 
electrons, Mn^+ in this material has four electrons in those five d-orbitals. Erom Hand’s rule, the spins are aligned 
ferromagnetically, and there are thus two possibilities for eg orbitals with either dx 2 -y 2 or d^z^-r^ being occupied. 
This represents the orbital degree of freedom in this Mott insulator, which can be modeled as pseudo-spins Tx^yy. 
The model Hamiltonian is 


= ^J4^T„(r)T^(r'), 


which is typically not SU ( 2 ) symmetric. With a long range orbital order, spin magnetism would be strongly affected 
by so called Jahn-Teller effect ( |Jahn and Teller 1937). 

Most p-orbital solid state materials, for example the semiconducting silicon and graphene, are actually weakly 
correlated. However, recent studies in one oxide heterostructure LAO/STO have found that correlated physics such 
as ferromagnetism emerges from the effective p-orbitals, where Px and Py are mimicked by dxz and dyz orbitals (the 
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FIG. 17 Interference measurement of inter-orbital coherence in the px + ipy superfluid (Kock et al. 
experimental protocol to prepare two copies of lattice condensates (red and blue). 




(a) shows the 


(b) shows the momentum distribution for 
the Px + ipy superfluid, (c) shows the atomic spatial distribution after ballistic expansion of the two condensates. The four 
Bragg peaks are labeled by 1-4. (d) shows the experimental observation of the interference pattern of the four Bragg peaks. 
The interference structure is along the z direction. 


degeneracy with dxy orbital is broken due to lack of out-of-plane inversion symmetry at the interface). In d-orbital 
systems, correlated physics usually emerges due to large Hubbard U interaction because of the tight confinement of 
these orbitals. The emergence of correlated physics in p-orbital systems on the other hand could be attributed to a 

2014b). In one dimension 


2013 


Li et al 


different origin, which is the quasi-one dimensionality (Chen and Balents 
at low filling, the magnetic susceptibility diverges as xid ^ where p is the occupation number per site. Even for 
infinitesimal interaction 7/, there is a strong interaction effect: the ratio of the interacting to free fermion susceptibility 
diverges, Xid/Xff ^ oo for p ^ 0. A general result for the free energy (per site) versus magnetization at low density 
is obtained to be 


F = 2pJeffFi 


M ksT 
P ’ Jefi 


-JhM^ 


( 120 ) 


where M is the magnetization (per site), Jh is the Hand’s rule coupling, Jeff is the effective antiferromagnetic coupling, 
and Fi[m,t] is the free energy per site of the one-dimensional antiferromagnetic chain, with reduced magnetization 
m and temperature t (this is kno wn from thermodynamic Bethe ansatz). The effective coupling Jeff is reasonably 
conjectured to scale as Jeff oc (Chen and Balents 2013). From the free energy, the Hand’s energy is dominant 


and favors a ferromagnetic state with sufficiently low density for arbitrarily weak Hand’s coupling Jh- A rigorous 


work (Li et al 2014b) studies the higher filling regime (but assumes no double occupancy), where a ferromagnetic 
ground state for p-orbital fermions is proved based on transitivity and non-positivity of the many-body Hamiltonian. 
Further studies are required to find out the boundary of the ferromagnetism in p-orbital fermions. 
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ZX 



t 2 g orbitals 



FIG. 18 Five (i-orbitals. In the presence of crystal field, the orbital degeneracy splits into two groups, eg and t 2 g- 


B. Synthetic orbital matter and material design 

In material science, design of materials for applications is an important subject. Recent developments involve 
engineering heterostructures with hybrid materials. For example oxide heterostructures such as LaAlOs/SrTiOs 
and GdTiOs/SrTiOs have been created and extensively studied. While the properties of many materials can be 
calculated within the density functional theory (DFT), this approach fails for ones of strong correlation for which 
d-orbital electrons typically play an important role. At the same time, these strongly correlated materials could have 
fascinating properties including important applications. High Tc superconductivity belongs to this class. Lots of 
efforts have been made in searching for materials with higher Tc, but there is no real improvement in the last two 
decades. Lack of reliable tools in predicting Tc leaves the design of high Tc superconductivity essentially to empirical 
trials, which are costly in both time and materials. Developing new tools to simulate strongly correlated materials by 
incorporating correlation effects in DFT has triggered tremendous interest but appears to be very challenging. 

To address the challenge of simulating correlated d-orbital electrons in classical computers, one alternative way 
is to create synthetic orbital matter with optical lattices and take it as a quantum orbital simulator. With this 
optical-lattice-based quantum orbital simulator, the ultimate procedure for material design would be—(1) conceive a 
particular design of materials; (2) determine the orbital configuration of the imagined material by quantum chemistry; 
and (3) apply cold atoms in optical lattices to simulate the properties. In such a way, we could explore the imagined 
quantum materials for desired properties, bypassing the often tedious chemical process of really fabricating them from 
electronic compounds. This would significantly speed up the material design and should help improve key quantities 
of great interest, for instance, the value of critical temperature Tc of superconductivity in future. Although the optical 
lattice experiment is still at a very early stage, with future developments, synthetic orbital matter in optical lattices 
could be extremely helpful to the design of real materials. 

Finally, we would like to point out that orbital degrees of freedom are found to play an important role for a 
vast majority of intriguing electronic quantum materials that condensed matter physicists have found since 1970s. 
Magnetic materials of spin only are an important class of systems that have been studied with great progress and 
remain to pose new challenges, such as frustrated magnets possibly showing spin liquid phases. In fact, the spin-only 
systems represent a small fraction of the world of real materials. Furthermore, past theoretical studies predicted exotic 
phenomena for model systems that have no spin but only orbital degrees of freedom. Such hypothetical models, which 
previously might have seemed too special and excessive, now become readily realizable with optical lattices. On this 
regard, using higher orbital bands of the optical lattice appears to open up a new front to explore orbital physics, 
both for understanding the electronic systems and for exploring artificial quantum orbital-only models that have no 
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prior analogue in solids. 


C. Many-body dynamics of high orbital atoms 


Coherent dynamics across different bands has been observed in many experiments ( 

Anderlini et al. 2007 Cheinet 

et al. 2008 Hu et al. 

2015 

Jona-Lasinio et al. |2003 

Sebby-Strabley et al. 2006 

Trotzky et al. 

2008 Zhdii et al. 2013). 

In particular the recent experiments ( 

Hu et al.\ 

2015 

Zhai et al. 

2013) have demonstrated fast coherent controllability 


of orbital degrees of freedom. These experimental developments open up possibilities of studying many-body dynamics 
of high orbital, where the observed Rabi-like oscillations between different bands can be affected by interaction. One 
particular example would be orbital Josephson effect, which has been studied for double-well potentials (Garcia-March 


et al. 2012 Garcia-March and Garr 2015 Garcfa-March et al 2011 Gillet et al. 2014). This effect has also been 


seen in numerical simulations of a dynamical procedure, proposed to detect the p ip BEG (Gai et al. 2012a Li 


et al. 2014a). 


The orbital Josephson effect is expected to be generic for various experimental setups for high orbital atoms. Here 


we consider the specific setup proposed to probe the complex order (see Sec. IV.E). Assuming all atoms condense, 
the dynamics is then approximately captured by a two-mode Hamiltonian, 


H = + h.c. 

T 9l ^Ki ^Ki + Ki —K2^ + ^2 ^Ki^K;i^K2^^2 

+ 9Z + h.c)j , 


( 121 ) 


where 5 k 1,2 condensed modes and the last term is a Umklapp process. Hollowing the treatment of 


Josephson effect developed for double-well Bose-Einstein condensates (Smerzi et al. 1997 Zapata et al. 1998), the 
dynamical state could be approximated by 


l^'W) = ;^ +V’2(i)&L) |0)- 


( 122 ) 


The corresponding time-dependent Gross-Pitaevskii equation is (Gai et al. 2012a) 

idt1pl{t) = Xlp2{t) + + g2\lp2?)lpl +253V’iV’2 > 

idttp2{t) = Xxpi{t) + {2gi\tp2\'^ + g2\'ipi\^)'ip2 + ■ 


(123) 


To make the dynamics more physical, one can rewrite the wavefunctions i^j{t) in terms of densities and phases as 

'>p2 

The equation of motion is most easily derived by constructing the Lagrangian, which takes the form, 

L = -pidtOi - P2dt02 - { 2 Av'/ 9 i /92 COs( 6»2 - 6 »i) + 253 / 01/02 cos( 2 ( 6*2 - 6 * 1 )) + giipl + pI) + g 2 plp 2 } ■ 

Erom Euler-Lagrangian equations, 

ddi ’ 


dL 

dpj 


= 0, 
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one gets 


dtpi = -dtp 2 = 2 \^piP 2 sin( 6>2 - 6 >i) 

+ 4^3Pip2 sin(2(6>2 - 6>i)), 

dtOi = -A. /^cos(6>2 - Oi) - 2gsp2 cos(2(6>2 - 6>i)) 

V Pi 

— 2^ipi — g 2 p 2 , 

dt02 = -A./^cos(Pi - O 2 ) - 2gspi cos(2(Pi - P 2 )), 

V P 2 

— 2^iP2 — 5 ' 2 Pi- 


To make a direct connection to Josephson effects, the number imbalance and phase difference are defined to be 
z = pi — p 2 and (j) = 0 i — 62 ^ whose dynamical evolution is governed by 


dtz = 2 ^A\/l — z‘^ sin(0) + g^{l — z^) sin 


2A^ 


dt4> = (2^1 - g2)z - ^2 cos(0) - 2^3 ^cos(20). 


(124) 


Compared with Josephson effects in double-well Bose-Einstein condensates (Smerzi et al. 1997 Zapata et al 


the key difference is that here we have sin(20) and cos(20)) terms which are generated by the Umklapp process gs. 
In the ground state, these terms give rise to the spontaneous time-reversal symmetry breaking. 

In the noninteracting limit, ^ 1 , 2,3 ^ 0, the Rabi-like oscillation with frequency 2A is easily recovered. In the linear 
regime, |z| ^ 1, the dynamics in z and 0 is simplified to 


dfZ ^ (2A + d^fs)J(/), 

dt6(l) ^ {2gi -g2-2X- 2gs)z, 


assuming 0 ^ 27r. This gives rise to oscillatory dynamics with a frequency 

CJreal = \/ (2A + 4^3)(2A — 2gi + 5^2 + 25 f 3 ), (125) 

which is the Josephson frequency for a real superposition state Px For the complex superposition px ± ipy^ 

expressing 0 in terms of fluctuation field 0 ^ | {5(j) <C 27r), the linear dynamics is 

dtz « -2gz{5(p- —), 

93 

dtS(p (2^1 - 52 + 253 ) 2 , 


which predicts a Josephson frequency 


^complex \/ 2^3 (2gi ^2 T 2 ^ 3 ), 


(126) 


with X/gs assumed to be small. In the Josephson effects, the frequency is different from that in non-interacting 
Rabi oscillations. This frequency difference is also seen in the numerical simulations based on Gross-Pitaevskii 


equations (Cai et al. 2012a) and Gutzwiller methods (Li et al. 2014a). 


The nonlinear effects of dynamics in Eq. (124) are expected to be more interesting, because of the sin(20) term, 
than the usual Josephson physics of double-wells. For example the analogy of self-trapping effect in double-wells would 
certainly exist in this orbital setting, and very likely would lead to new possibilities beyond the standard double-well 
Josephson effect. Details of such orbital Josephson effects call for further theoretical and experimental investigations. 


D. Relation to spin-orbit coupled quantum gases 


Orbital degree of freedom can certainly be mapped to pseudo-spins. In doing so, spin-orbit couplings of certain 


types usually arise naturally due to the spatial nature of orbitals (Belemuk et al. 2014 Li et al. 2013 Liu et al. 2014 


2010 Sun et al. 2012a|b Zhou et al. 2015). The tunneling Hamiltonian of orbital models mixes different orbitals. In 


particular mixing of different parities could lead to non-trivial effective spin orbit couplings and consequent topological 
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properties. Mixing of s and p-orbitals in a ladder system ( 

Li et al. ^ 

^013) closely mimics the one dimensional spin 

orbital coupling recently engineered in cold gases by Raman transitions 

(Galitski and Spielman 2013 

Lin et al. 

20111 

Such sp orbital mixing is recently achieved in a shaken lattice experiment of ^^^Cs Bose-Einstein condensates (Parker 


et al. 2013) and a similar band structure with double minima like the spin-orbit coupled case is indeed obtained. 


Mixing of p and d-orbitals gives rise to the phases of topological semimetal and topological insulator (Sun et al. 2012b). 
One recent work shows that mixing of p-orbitals in spin imbalanced fermions leads to topological superconductivity 
with novel features (Liu et al, 2014). We note however that some of these novel predictions are made for fermionic 
species of atoms, whereas the high band experiments have been explored only for bosons so far as to this time. Further 
experimental developments are expected. 

With strong repulsion, particles could form Mott states with the charge degrees of freedom frozen. The orbital 
ordering in Mott states is then described by super-exchange interactions of orbitals, which typically depend on the 
orientation of links. T his orientation dependent orbital super-exchange gives rise to novel pseudo-spin models such as 
quantum 120^ model (Zhao and Liu 2008 Wu 2008b) (see Eq. ([^). 

With spin-orbit couplings, many interesting quantum phases such as skyrmions and topological states have been 
investigated. The connection of orbital physics to spin-orbit coupling suggests possibilities of novel orbital states. 
One reason to study spin-orbit coupled physics in orbital systems (with atoms loaded into higher bands) is that there 
appears no additional heating in this system, in contrast with the heating challenge faced by engineered spin-orbital 
couplings by the advanced Raman laser technique. In this regard, orbital physics provides an alternative platform to 
investigate spin-orbit coupled phenomena, which is a direction worth future exploration. 


E. Periodic driving induced orbital couplings 


In recent optical lattice experiments 

(Aidelsburger et al. 

2013 

Jotzu et al. 2014 Miyake et al. 

2013 

Niu et al. 

2015 Parker et al. 2013 

Struck et al. ' 

2013 Weinberg et al. 

2015 

), periodically driven systems have been developed 


with a motivation to create exotic atomic phases. In such systems time reversal symmetry is explicitly broken. With 
the driving frequency matching band gaps, energetically separated orbital bands can be efficiently coupled. 

Here we use one example to demonstrate the key idea of using lattice shaking to induce/control orbital couplings. 
Consider a one dimensional shaking lattice as implemented in experiments (Parker et al. ] |M^ . The time-dependent 
optical potential of this lattice reads 


V{x, t) = Vb cos [k{x — xo{t))] , (127) 

with xo{t) a periodic function, xo{t) = Xq sin(27rt/T). Taking Xq = 0, we have a static lattice potential where s 
and p orbital bands are decoupled and well separated by an energy gap. With weak driving, we have V{x,t) « 
Vb [cos{kx) kxo{t) sm{kx)]. The time-dependent term introduces an effective coupling between s and p orbitals, 
approximately given by 

Xsp = kVoXo{t) j dx sin(/cx)rc*(x)rcp(x), (128) 

with Wu{x) the orbital wavefunction. With frequency 27r/T matching the band gap, the system is approximately 
described by a static two-band model with s and p orbitals coupled, under a rotating wave approximation. 

It appears natural to engineer orbital couplings by lattice modulation/shaking techniques. But the problem is that 
heating effects are fundamentally unavoidable in periodically driven quantum systems. Since periodic driving breaks 
time translational symmetry, energy is no longer a conserved quantity. It follows that driven systems (assuming 
ergodicity) at long time would necessarily be described by infinite temperature ensemble. Nonetheless, there could be 
long lifetime transient states that manifest interesting topological features. This requires more careful treatment of 
quantum dynamics than just solving for the ground states of effective static Hamiltonians. One way out is to combine 
with dissipation. Driven-dissipative orbital models may exhibit steady quantum many-body states with interesting 
topological properties. This is worth future exploration. 


F. Open questions 

For bosons, firstly, it remains open how to experimentally reach the Mott insulator phases of the p-band and 
study the p-band superfluid-Mott insulator transition. The current experiments at Hamburg are performed with a 
two-dimensional checkerboard lattice and a relatively shallow harmonic trap in the third dimension. Introducing an 
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additional optical lattice potential in the third dimension is required to access the Mott regime. Unfortunately that 
would also increase the on-site interaction between p-orbital bosons, which leads to faster decay (Hemmerich 2014). 

Secondly, it is intriguing to find out what type of new topological defects, other than vortices, may possibly occur 
in the staggered Px zb ip^-orbital Bose-Einstein condensate. The state breaks not only U(l) but also other interesting 
symmetries that are usually not broken in other conventional Bose condensates, including for example, time-reversal, 
lattice translational and rotational symmetries. On the general ground of broken symmetries, new classification of 
topological defects is expected but remains unknown. 

For fermions, the stability of the p and higher orbital bands is protected by Fermi statistics, if the experimental 


(Kock et al. 

2015 

Muller et al. 

2007 


Olschlager et al. 


2013 


2012 


Wirth et al. 


2011). Nevertheless, this 


approach would require a higher density of fermions, which in turn requires a higher efficiency of cooling fermions 
down to degeneracy. The recent breakthrough in the Rice experiment of fermions on lattice (Hart et al. 
promising for studying the higher orbital bands. 


2015) is 
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Appendix A: Tree level estimate of couplings in effective field theory for p-orbital bosons 


In this appendix, the coupling constants in the effective field theory (Eq. ( [3Q| )) are related to a microscopic model. 
We start with the contact interaction for a 3D Bose gas, which reads 

271 a. 


Unt,= 




(Al) 


where '^(x) is a bosonic field operator, m is the mass of atoms and is the 3D scattering length. With bosons loaded 
into the p-band of a 2D lattice that has band minima at = (tt, 0) and = (0, tt), the field operator is expanded 
by the low energy modes as (Li et al. 2014a) 


1p{x) = J 




{^Qx+qe 




MQ,+q(x) 




(A2) 


where 5q +q is the annihilation operator for a Bloch mode near the band minimum and uq +q(x) is the cor¬ 
responding periodic Bloch wavefunction. At tree level, the high energy modes are integrated out and the resulting 
renormalization of the low energy theory is neglected. Then the interaction is written in terms of these low energy 
modes as 




27 Tash^ 

m 


[nbi 


(27r)2 


g-i(qi+q3-q2-q4)-x 


{^Qx+qjQ-+‘J2^Qx+q3^Q»+‘i4“Qx+qiW“Qx+q2W«Qx+q3W“Qx+q4W + ^ Qy 

&L+qi^Qx+q2^Q,+q3^Q«+q4^Q=.+qi(x)^iQ.+q3(x)WQ^+q3(x)MQ^+q,(x) 

^Q.+qjQ«+‘l2^L+q3^Q.+'l4“Q.+qlx)«Q,+q3(x)WQ,+q3(x)«Q„+q,(x) + O Qj,} . 


(A3) 
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We can rewrite x = R + x', where R is the position vector of lattice sites and x' centers over one unit cell. With 

— q^) (a jg the lattice constant), we get 




-I 


fi (2^)^ 


(27r)^5(qi + Qs - q2 - q4) < Y1 ^2, q3)^Q„+qi^Q„+q2^Q^+q3^Q^+q4 


a,P=x,y 


where 


<73(qi,q2,q3) [^Q,,+qi^Q,+q3^Q«+q2^Q«+q4 + Qa. ^ Qy]} > 


5a;a;(qi,q2,q3) = J C^^x'uQ^+q^(x')UQ^+q3(x')UQ^+q3(x')UQ^ 


ma 

9yy ^2, = J (x^UQ^+q^ (x^U^j^+q^ (x')UQ, 

Wy(qi,q2>q3) = flya:(q3>qi + q3 - q2>qi) 

y t^^x'u^^+q^(x')UQ^+q^(x')u^^+q3(x')UQ^+q3+q3-q3(x'), 


+qi+q3-q2(^ )’ 
+qi+q3-q2(^ )’ 


5'3(qi,q2,q3) = 


ma^ 

ATTa.h? 


ma^ 


J <^^x'uQ^+q^(x')UQ^+q3(x')UQ y+q2 +qi+q3-q2(^ )• 

Neglecting the momentum dependence of gai 3 and ^ 3 , the derived couplings simplify to 

271 a. 


9xx 9yy 


9xy 9yx 


ma^ 

Airagli? 

ma^ 


j dVluQ^x')!"^, 


271 a. 


93 = 


ma^ 


J C^^x'uQJx')*^UQ^(x')^ 


(A4) 


(A5) 


The calculation of and Kj_ is straightforward at tree level, and they are estimated to be K\\ = — ^ ^2 ^p(k) |k^Q^, 


and = — ^^£’p(k)|k^Q^, with £’p(k) the dispersion of the p-band. 
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